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1 .0  Objectives 


The  objective  of  the  effort  as  described  in  the  proposal  to  DARPA/DSO  was  to  identify 
and  assess  the  primary  technical  barriers  and  risks  that  must  be  overcome  to  establish 
fast  computational  methods  for  uncertainty  analysis  in  large  models  of  interconnected 
components  where  the  shape  of  system-level  dynamics  is  essential.  This  objective  was 
addressed  through  the  following: 

1 .  Validate  the  technical  hypothesis  that  measure-theoretic  algorithms  (for 
network  congestion  in  large  systems  of  interconnected  model  components) 
could  be  used  to  radically  accelerate  the  numerical  simulation,  identification 
of  dynamic  stability  boundaries,  model  reduction  and  propagation  of 
uncertainty  during  product  design. 

2.  Apply  preliminary  results  to  selected  challenge  problem  applications  in  (i)  the 
design  of  organic  molecules  where  shaping  the  conformational  dynamics  is 
essential  to  engineering  their  exploitable  properties;  and  (ii)  the  control  of 
afterburner  combustion  instabilities  in  gas  turbine  propulsion  systems 
designed  for  advanced  military  aircraft. 

3.  Champion  advanced  computational  methods  for  the  management  of 
uncertainty  in  large,  interconnected  systems  where  the  control  of  dynamics  is 
essential.  Engage  international  communities  of  experts  in  mathematics  and 
engineering  to  sharpen  and  articulate  an  actionable,  mathematically  rigorous 
approach  for  the  development  of  Analytic  Systems  Engineering  with  the 
potential  of  broad  impact  on  both  the  military  and  commercial  sectors. 

Table  1  describes  our  results  and  accomplishments  during  the  contract  period  on  each  of 
the  required  contract  task  elements  designed  to  achieve  these  objectives. 
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Table  1  -  Task  elements  for  contract  F49620-03-C-0035 


Task 

Description 

Status 

Task  A 

Establish  the  correct  technical  framework  for  model- 
based  uncertainty  analysis  using  the  interconnection 
(graph)  structure  between  components. 

Completed 

(see  Publications 

5,  6,  and  7) 

Task  B 

The  analytical  framework  described  in  Task  A  will  be 
used  to  access  (i)  the  potential  of  the  graph-theoretical 
methods  for  accelerating  computational  methods  for 
estimating  the  propagation  of  uncertainty  measure  in 
interconnected  model  systems;  and  (ii)  the  technical 
obstacles  and  risks  involved  in  designing, 
implementing,  and  testing  numerical  algorithms  for 
model-based  uncertainty  management  in  the  design  of 
complex  systems. 

Completed 

(i)  Addressed  in 
Publications 

5,6,7 

(ii)  Addressed  in 
Publications  3,4 

Task  C 

In  collaboration  with  the  DARPA/DSO  program 
manager,  specific  model  systems  will  be  selected  to 
investigate  where  the  connection  structure  arises  from 
the  physical  proximity  between  the  subunits.  The 
primary  application  domain  will  be  the  design  and 
control  of  conformation  dynamics  in  large  organic 
molecules,  although  related  problems  will  be 
considered. 

Addressed  in 
Section  3.7.1 
and  Publication  8 

Task  D 

In  collaboration  with  the  DARPA/DSO  program 
manager,  specific  model  systems  will  be  selected  to 
investigate  where  the  connection  structure  arises  from 
the  virtual  relationships  between  the  subunits.  This 
will  be  a  model  of  combustion  instability  that  occurs  in 
the  augmentor  in  a  high  performance  gas  turbine 
engine  for  advanced  military  aircraft. 

Addressed  in 
Publication  2 

Task  E 

For  the  selected  model  systems  selected  in  Tasks  C-D 
above,  we  will  establish  the  computational  complexity 
bounds  as  they  scale  with  model  size  (e.g.,  as  size  of  a 
model  protein  length  or  with  the  number  of  coupled 
interacting  modes  in  a  combustion  system). 

Addressed  in 
part  in  Monte 
Carlo, 
Polynomial 
Chaos,  and  DNA 
work 

Task  F 

Host  a  DARPA/DSO  sponsored  Workshop  involving 
recognized  experts  in  various  areas  of  Mathematics  and 
Engineering  to  evaluate  the  state-of-the-art  and 
potential  for  fast  techniques  supporting  uncertainty 
analysis  in  large,  multi-scale  dynamical  systems  with 
applications  in  computational  biology,  design  of 
advanced  weapon  systems  and  information  technology. 

Completed 

Workshop  held 
Jan. 15-16,  2004 
atUTRC 
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2.0  Work  Performed 


A  multi-faceted  effort  into  the  propagation  of  uncertainty  in  complex  networked  systems 
was  conducted.  This  effort  included: 

1 .  Graph  theoretic  decomposition  schemes  for  converting  a  large  monolithic 
system  of  differential-algebraic  equations  into  a  collection  of  weakly  coupled 
subsystems  of  strongly  coupled  equations.  Computational  benefits  are 
achieved  by  solving  less  stiff  systems  in  parallel. 

2.  Spectral  balance  methods  for  large  coupled  limit  cycling  systems  for 
applications  in  aeroengine  augmentor  instability  control. 

3.  Determination  of  density  mapping  techniques  for  propagation  probability 
density  functions  through  large  systems  of  interconnected  component  models. 

4.  Uncertainty  propagation  schemes  utilizing  polynomial  chaos  methods  for 
nonlinear  systems  with  nontrivial  dynamics. 

5.  Symmetry-breaking  schemes,  where  we  deliberately  introduce  spatial 
variations  in  the  system  parameters  in  order  to  change  the  stability  properties. 

6.  Uncertainty  propagation  in  dynamical  systems  involved  in  computational 
chemistry  and  molecular  modeling. 

7.  DNA  modeling,  where  we  investigated  dynamical  systems  where  the 
connection  structure  arises  from  the  physical  proximity  between  the  subunits. 

8.  Uncertainty  in  the  dynamics  of  conservative  maps,  focusing  on  the  standard 
map  and  a  discrete  Duffing  oscillator  in  R2. 


3.0  Accomplishments  and  New  Findings 

In  the  paragraphs  that  follow,  we  detail  the  key  technical  work  done  in  conjunction 
with  this  contract.  The  work  touched  on  a  number  of  seemingly  diverse  topics,  that, 
when  taken  in  their  entirety,  form  the  basis  for  an  analytic  understanding  of  uncertainty 
propagation  through  large,  interconnected  systems. 

3.1  Literature  Search 

An  extensive  literature  search  was  performed  in  the  area  of  uncertainty  propagation 
with  emphasis  on  dynamical  systems  that  can  be  modeled  as  interconnected  structures. 
Information  was  gathered  in  the  areas  of  large,  multiscale,  highly  interconnected  systems 
where  dynamics  plays  an  important  role.  Papers  were  identified  covering  graph  theoretic 
methods  to  preserve  structure  dictated  by  the  local  dynamics.  Papers  on  spectral  methods 
for  invariant  measure  computations  in  dynamical  systems  were  collected.  Monte  Carlo 
methods  were  reviewed.  Investigations  of  newer  techniques  such  as  polynomial  chaos 
and  stochastic  finite  element  analyses  were  collected.  Existing  methods  in  molecular 
modeling  were  reviewed.  Literature  on  methods  of  harmonic  analysis  and  describing 
functions  in  control  and  dynamical  systems  were  reviewed.  A  study  of  existing  software 
products  that  purport  capabilities  in  uncertainty  propagation  was  done.  Most  such 
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products  are  in  the  area  of  reliability  analysis  and/or  probabilistic  optimization.  A 
repository  of  all  identified  documents  is  being  maintained  and  all  documents  are  available 
to  all  members  of  the  technical  and  programmatic  team. 


3.2  Uncertainty  Propagation 

Uncertainty  analysis  is  a  research  topic  that  has  received  much  attention  in  recent 
years.  Increased  use  of  physics-based  models  in  the  study  of  dynamical  systems,  in  a 
wide  range  of  applications,  calls  for  model  prediction  analyses  and  quantification  in 
terms  of  uncertainties  in  descriptions  and  operating  environments.  In  this  section  we 
focus  on  several  advanced  methods  for  propagating  and  evaluating  the  effects  of 
uncertainty  in  complex  dynamical  model  parameters  and  their  initial  conditions. 

Polynomial  Chaos  (PC)  (i.e.  Stochastic  Finite  Elements)  is  an  analytical  approach 
based  on  expansions  of  uncertain  quantities  in  terms  of  prescribed  random  basis 
functions.  This  methodology  has  received  significant  attention  during  recent  years  and 
has  been  applied  to  uncertainty  propagation  in  complex  dynamical  systems,  which  are 
described  by  partial  differential  equation  models.  It  has  been  demonstrated  that  for  a 
certain  class  of  problems  arising  in  fluid  mechanics,  finite  elements  and/or  chemical 
systems,  PC  can  be  considerably  (up  to  several  orders  of  magnitude)  faster  than  Monte 
Carlo  (MC)  or  similar  methods.  Furthermore,  the  analytical  representation  in  the  PC 
method  can  be  of  great  benefit  in  system  analysis. 

3.2.1  Monte  Carlo  Methods 

In  the  framework  we  proposed  for  uncertainty  analysis  of  large  nonlinear  systems, 
the  most  important  objects  are  invariant  measures.  These  can  always  be  computed  using 
Monte  Carlo  methods,  usually  at  a  high  cost.  We  studied  application  of  Monte-Carlo 
methods  in  a  set  of  dynamical  systems  with  known  asymptotic  dynamics.  Uncertain 
inputs  and/or  parameters  were  randomly  sampled  against  some  predetermined  probability 
density  functions  and  the  model  was  run  many  times.  A  wide  variety  of  sampling 
techniques  were  investigated,  including  pure  random  selection,  randomized  quasi-Monte 
Carlo  techniques,  and  Latin  hypercube  methods.  Approximately  ten  different  probability 
density  functions  were  identified  to  sample  from  in  order  to  allow  for  a  generalized 
technique.  Additionally,  the  importance  sampling  methods  were  investigated,  as  well  as 
Markov  chain  Monte  Carlo  methods. 

It  is  well  known  that  Monte  Carlo  methods,  while  guaranteed  to  work  over  the  long 
run,  are  poorly  suited  to  high  dimensions  (large  numbers  of  uncertain  inputs  and/or 
parameters).  Additionally,  since  numerous  evaluations  of  the  overall  model  are  required, 
these  methods  suffer  large  penalties  if  the  model  evaluation  is  costly,  either  in  time  or  in 
dollars.  Finally,  there  is  no  appropriate  way  to  allow  the  inherent  structure  induced  on  the 
problem  by  the  physical  nature  of  the  problem  to  be  useful  in  reducing  the  computational 
complexity.  Monte  Carlo  techniques  are  the  epitome  of  “black  box”  methods,  were  no 
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understanding  of  the  model  is  used.  Monte  Carlo  methods  can  be  considered  for  lower 
dimensional  base-line  analyses,  to  measure  improvements  over  other  candidate 
approaches.  Additionally,  the  method  can  be  used  in  the  cases  when  polynomial  chaos 
expansion  fails  due  to  the  compact  nature  of  the  asymptotic  dynamics,  to  provide  a 
simulation  tool  for  measure  propagation  in  localized  dynamics. 


3.2.2  Polynomial  Chaos 

In  1938,  Wiener[l]  investigated  an  approach  he  called  “homogeneous  chaos”.  This 
work  was  expanded  by  Ghanem  and  Spanos[2].  The  idea  is  to  project  the  variables  of  the 
problem  onto  a  “stochastic  space”,  spanned  by  a  set  of  complete  orthogonal  polynomials 
that  are  functions  of  the  random  variables.  The  goal  is  to  eliminate  the  large  sample 
requirments  of  Monte  Carlo  methods,  thereby  dramatically  reducing  the  computational 
requirement  in  the  study  of  uncertainty  propagation. 

In  Wiener’s  original  work,  his  random  variables  were  assumed  normally 
distributed,  and  he  used  the  Hermite  polynomials  as  his  basis  functions.  Xiu  and 
Kamiadakis[3]  investigated  alternative  pairs  of  random  variables  and  orthogonal 
polynomial  basis  functions,  such  as  gamma  distributed  random  variables  with  Laguerre 
polynomials. 

Our  efforts  focused  on  developing  an  analytical  understanding  of  these  methods  and 
identifying  the  class  of  problems  where  they  may  show  promise  and  classes  of  problems 
where  polynomial  chaos  methods  are  not  suitable  for  uncertainty  propagation.  We  have 
studied  a  hierarchy  of  models  starting  from  simple  linear  models  with  parametric 
uncertainty,  to  limit  cycling  models  with  additive  noise. 

For  a  large  class  of  smooth  dynamical  systems,  these  polynomial  chaos  methods  are 
effective  in  markedly  reducing  the  computational  burden  inherent  in  the  Monte  Carlo 
methods  mentioned  above.  In  our  experience  we  have  found  two  to  three  orders  of 
magnitude  reductions  in  computational  effort.  However,  we  have  also  found  limitations. 

It  was  discovered,  following  earlier  work  by  Kamiadakis,  that  polynomial  chaos 
expansions  can  break  down  if  the  asymptotic  dynamics  resides  on  a  compact  attractor  for 
all  (compactly  supported)  uncertain  parameter  values.  A  solution  to  this  is  proposed  to  be 
the  development  of  spectral  element  methods  tuned  to  the  asymptotic  dynamics  of  the 
nonlinear  system  under  study.  In  addition,  multiscale  spectral  element  methods  could  be 
used  to  provide  multiscale  uncertainty  propagation  methods.  This  way,  the  substantial 
speed-up  achieved  by  polynomial  chaos  methods  can  be  preserved  for  asymptotically 
compact  uncertain  attractors.  In  addition,  the  above  described  utilization  of  the  graph 
structure  of  the  problem  is  used  to  render  polynomial  chaos  methods  useful  for  a  class  of 
large  networked  systems,  where  a  brute-force  expansion  might  fail  due  to  a  large  number 
of  uncertain  parameters. 

The  truncated  polynomial  chaos  expansion  has  a  structure  of  an  A-dimensional 
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dynamical  system,  where  N  is  mxn,  m  is  the  dimension  of  the  phase-space  for  the 
dynamical  system  under  study  and  n  is  the  order  of  truncation.  We  have  proven  some 
properties  of  the  resulting  dynamical  system  for  the  cases  discussed  above.  In  particular, 
it  can  be  shown  that  for  a  compact  uncertain  attractor,  the  polynomial  chaos  expansion 
asymptotically  diverges  due  to  the  fact  that  the  output  measure  has  its  density  in  a  “bad 
space”. 


References  for  this  section:  (Note:  Each  section  has  its  own  references,  renumbered, 
starting  at  1 .) 

[1]  Wiener,  N.,  “The  homogeneous  chaos”,  Amer.  J.  Math.,  60,  pp.  897-936,  1938 

[2]  Ghanem,  R.,  and  Spanos,  P.,  Stochastic  Finite  Elements A  Spectral  Abroach,  Springer- Verlag, 
New  York,  1991. 

[3]  Xiu,  D.  and  Kamiadakis,  G.E.,  “The  Wiener-Askey  Polynomial  Chaos  for  Stochastic  Differential 
Equations”,  SIAM  J.  Sci.  Comp.,  Vol.  24,  No.  2,  pp.  619-644,  2002. 

3.2.3  Density  Mapping 

Enhancing  engineering  performance  and  productivity  through  systematic  and 
integrated  designs  that  account  for  the  effects  of  uncertainty  is  a  key  economic 
priority[  1  ] .  A  unified  paradigm,  referred  to  as  "Analytical  Systems  Engineering",  recently 
formulated  at  the  United  Technologies  Research  Center,  has  identified  uncertainty 
propagation  through  networks  of  nonlinear  components  as  an  essential  component.  This 
view  is  also  echoed  in  [2]  or  in  [3]  (at  the  end  of  this  secton),  where  large-scale, 
interconnected  systems  are  designed  using  model-based  techniques  that  employ 
uncertainty  descriptions  explicitly.  Since  all  system  models  have  varying  levels  of 
uncertainty  [4],  designers  often  use  large  safety  margins,  which  result  in  more  complex 
and  expensive  systems  [5].  As  a  result,  a  natural  path  in  modem  systems  design  is  to 
make  decisions  on  the  best  system  structure  from  the  perspective  of  greater  robustness  to 
uncertainty. 

Before  moving  towards  design,  addressing  uncertainty  analysis  in  a  model-based 
complex  engineering  systems  framework  is  important.  This  has  traditionally  been 
addressed  by  Monte  Carlo  (MC)  like  methods.  This  classic  approach  (see  [6])  employs  a 
large  number  of  simulations  with  a  random  selection  of  variables  from  their  prescribed 
distribution  (parametric  or  empirical). 

Unfortunately,  uncertainty  propagation  techniques  using  MC  methods,  even  in  their 
advanced  forms,  do  not  scale  well  with  system  size.  The  natural  choice  in  such  cases  is  to 
break  the  large  system  into  pieces.  More  precisely,  the  physics-based  models,  which  are 
often  converted  to  monolithic  systems  of  uncertain  nonlinear  differential/algebraic 
equations,  are  to  be  decomposed  using  graph  decomposition  methods  into  subsystems 
evolving  on  different  time  scales,  as  mentioned  in  [7]. 
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We  find  that  the  time  scale  separation  can  be  exploited  to  increase  computational 
efficiency  when  propagating  input  uncertainty  in  a  subsystem-by-subsystem  manner.  This 
approach  has  been  also  advocated  in  [8],  where  arbitrary  interconnections  of 
multivariable  systems  (represented  either  in  a  continuous  or  discrete  form)  with  nonlinear 
or  linear  dynamics  (nonlinear  time  varying,  distributed  linear  time  invariant  or  lumped 
linear  time  invariant)  are  decomposed  into  aggregate,  strongly  connected  subsystems. 
Following  this  procedure  each  subsystem  is  addressed  using  "the  minimum  set 
technique"  and  transformed  into  a  typical  feedback  interconnection. 

Alternatively,  in  [9],  a  Recursive  Projection  Method  (RPM)  is  developed  to  solve 
nonlinear  parameter  problems  for  which  the  convergence  is  achieved  for  certain 
parameter  values.  The  correction  applied  for  divergent  domains  involves  Newton's 
integration  method.  RPM  provides  reliable  results,  when  the  number  of  divergent  modes 
is  small  comparative  to  the  system's  dimension.  The  algorithm  has  been  successfully 
demonstrated  on  folds,  bifurcations  and  unstable  system  branches.  This  approach  is 
believed  by  its  authors  to  greatly  accelerate  iteration  convergence. 

Similarly,  extensions  from  the  backward  Euler  formula,  conventionally  used  to 
obtain  a  system  of  nonlinear  algebraic  equations  from  an  original  system  of  nonlinear 
algebraic  differential  equations,  can  be  grouped  under  the  waveform  relaxation  (WR) 
method.  This  technique  has  been  successfully  employed  in  [10]  to  address  large  scale 
systems,  such  as  integrated  circuits.  The  iterative  WR  method  decomposes  the  system  in 
several  dynamical  subsystems,  which  are  independently  analyzed  for  the  entire  time 
interval.  This  method  comes  with  sufficient  convergence  guarantees,  also  revealed  in 
[10]- 


In  [1 1]  graph  theory  is  also  used  in  the  context  of  autocatalytic  networks/sets  to 
classify  the  uncertainty  of  the  network  and  predict  its  influence  over  short  and  medium 
time-scales.  Therefore,  looking  at  how  networks  evolve  with  time,  more  precisely 
looking  at  their  dynamics,  it  can  be  beneficial  from  a  computation  speed  perspective.  It  is 
essential,  when  dealing  with  irreducible  (i.e.  strongly  connected)  graphs  (e.g.  dynamical 
subsystems),  to  perform  related  computations  in  corresponding  time  scales  [12]. 

Finding  a  subset  of  the  state  space  of  a  dynamical  system  where  typical  trajectories 
stay  longer  before  entering  different  regions,  conventionally  called  almost  invariant  sets 
[13],  entitles  a  macroscopic  behavior  analysis  of  dynamical  systems  networks.  In  effect 
this  separation  permits  individual  sub-system  uncertainty  propagation  studies.  Algorithms 
for  containment  of  such  almost  invariant  sets  are  presented  in  [13]. 

To  address  slow  convergence  rate  and  clustering  issues  associated  with  Monte 
Carlo  methods,  new  methods  such  as:  polynomial  chaos  [14]  (i.e.  stochastic  finite 
elements),  stochastic  surface  response  methods  [15]  and  probabilistic  collocation 
methods  [16]  were  used  with  significant  success.  Complementing  these  approaches,  the 
propagation  of  uncertainty  in  the  distribution  of  initial  conditions  for  a  network  of 
dynamical  systems  can  be  studied  in  an  exact  manner  using  Liouville's  equation  as  in 
[17]. 
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A  bottleneck  in  the  application  of  uncertainty  propagation  methods  in  a  subsystem- 
by-subsystem  way  is  associated  with  the  existence  of  correlated  signals  due  to  typical 
parallel  or  feedback  connections  found  in  conventional  subsystems.  To  address  this  issue, 
the  propagation  of  uncertain  inputs  through  interconnections  of  dynamical  systems,  with 
simple  asymptotic  behavior,  was  studied.  The  method  of  choice  is  a  discrete  density 
mapping,  analogous  to  the  input-output  Perron-Frobenius  operator. 

The  advocated  method  assumes  that  a  large  system  is  broken  into  components, 
which  evolve  on  different  timescales  or  have  simple  attractor  structure.  Then,  uncertainty 
propagation  methods  are  used  to  map  the  input  distribution  through  components  and 
arrive  at  a  stationary  density  for  the  states/outputs  of  interest.  In  a  chain  topology  only  the 
need  for  a  transfer  operator  is  observed.  This  is  in  contrast  with  the  parallel  connection, 
which  adds  complexity  through  the  necessity  of  summing  correlated  signals.  In  the  case 
of  feedback  structures  this  complexity  is  further  augmented  by  the  requirement  for 
convergence  of  the  iterations. 

An  electronic  version  of  this  work  in  its  full  depth  is  included  here.  The  same  paper 
can  be  found  in  the  appendix  of  the  printed  report  and  is  entitled:  Propagation  of 
Uncertain  Inputs  Through  Networks  of  Nonlinear  Components. 
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3.3  Graph  Decomposition 

Complex  physical  systems  often  possess  an  inherent  structure  from  which  they 
derive  the  robustness  of  their  behavior.  For  example,  the  DNA  molecule  has  a  strong 
backbone  structure  that  allows  global  modes  [1]  and  gene  regulatory  networks  have  a 
hierarchical  structure  that  allows  complex  gene  expressions  [2].  Similarly,  in  complex 
engineering  systems  such  as  coordinated  group  of  vehicles,  sensor  and  communication 
networks  are  designed  by  interconnecting  subsystems  in  a  systematic  manner.  There  is 
significant  interest  in  robust  design  methods  for  complex,  interconnected  dynamical 
systems.  The  design  of  such  systems  can  be  significantly  accelerated  by  using  models  for 
robustness  assessment  and  redesign.  This  requires  understanding  how  uncertainty  in 
various  parameters  and  the  system  model  itself  affects  the  model  outputs.  Computational 
methods  such  as  Monte-Carlo,  polynomial  chaos  and  various  modifications  of  these  can 
be  used  to  propagate  probabilistic  information  from  the  inputs/parameters  to  the  outputs 
of  a  static  or  dynamical  system.  Unfortunately,  the  computational  effort  for  these 
methods  scales  poorly  with  system  complexity  and  there  is  a  strong  need  for  efficient 
alternatives.  A  promising  direction  is  the  decomposition  of  a  complex  system  using  graph 
theoretic  methods  followed  by  block-by-block  propagation  with  iterations  when 
necessary. 

Graph  theoretic  methods  are  used  widely  in  control  theory  [3],  computer  science, 
and  network  systems.  More  recently,  graph  theoretic  techniques  have  been  employed  in 
conjunction  with  set  oriented  numerical  methods  for  computing  the  almost  invariant  sets 
of  dynamical  systems  [4].  There  is  a  close  connection  between  Markov  chain  theory  and 
graph  theory  [5].  We  refer  the  reader  to  standard  text  books  and/or  review  papers  such  as 
[6],  [7]  for  basic  notions  from  algebraic  graph  theory  and  graph  algorithms  such  as 
spanning  and  induced  subgraphs,  cycles,  connectedness  of  undirected  graphs,  strongly 
connectedness  and  topological  sorting  of  directed  graphs,  adjacency,  incidence  and 
Laplacian  matrices,  the  significance  of  the  eigenvalues  of  these  matrices,  reducibility  and 
irreducibility  of  the  adjacency  matrix,  depth  first  search  and  breadth  first  search. 

The  aim  of  this  effort  was  to  introduce  graph  theoretic  methods  for  decomposing  a 
complex  dynamical  system  into  hierarchically  and/or  weakly  connected  subsystems  for 
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the  purpose  of  uncertainty  propagation  in  systems.  The  structure  resulting  from  graph 
decomposition  with  series,  parallel  or  feedback  connections  among  subsystems  can  be 
exploited  for  efficient  methods  for  computation  of  uncertainty  propagation  [8].  In  the 
paper  included  here  (electronic  version,  and  in  the  appendix  in  the  print  version  of  this 
report),  entitled:  Graph  Decomposition  Methods  for  Uncertainty  Propagation  in 
Complex,  Nonlinear  Interconnected  Dynamical  Systems,  we  briefly  review  a  structural 
decomposition  algorithm  [3],  [9]  which  employs  the  so-called  equation  graph  of  a 
dynamical  system  and  identifies  the  subsystems  and  the  associated  hierarchy. 
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3.4  Spectral  Balance 

Many  industrial  flows  involve  complex  interactions  of  acoustic  waves,  vorticity, 
fuel  transport,  and  chemical  reactions.  The  control  objective  often  is  to  create  beneficial 
non-equilibrium  dynamics  with  control.  Examples  include  control  of  flow  separation  and 
mixing  enhancement.  Our  effort  here  focused  on  a  frequency  domain  framework  for 
analysis  and  non-equilibrium  control  design  for  a  large  class  of  models  of  physical 
phenomena  involving  multiple  oscillatory  modes  coupled  through  nonlinear  terms  and/or 
transport  delay,  that  are  driven  by  broad-band  disturbances.  While  motivated  by  specific 
problems  arising  in  military  aeroengines,  the  methods  developed  are  applicable  to  a  large 
class  of  distributed  dynamical  systems  involving  oscillatory  dynamics  with  nonlinear 
cross-coupling,  saturated  nonlinearities,  transport  delay,  and  broad-band  disturbances. 
The  spectral  balance  framework  that  we  developed  generalizes  the  standard  harmonic 
balance  and  Gaussian  signal  balance  in  feedback  systems  [1],  [2],  The  framework  was 
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introduced  and  illustrated  in  an  example  of  a  nonlinear  system  that  exhibits  noise-induced 
transitions  between  two  stable  equilibria.  The  example  presented  was  a  scalar  model  with 
cubic  nonlinearity  after  a  pitchfork  bifurcation  driven  by  a  broad-band  disturbance.  An 
approximate  and  iterative  spectral  balance  of  the  constant  and  broad-band  signals 
(including  determination  of  equilibria)  was  solved.  The  solution  of  this  approximate 
spectral  balance  was  used  to  reformulate  the  original  model  using  a  loop  transformation 
so  that  an  iterative  procedure  for  finding  the  spectrum  of  the  output  converges  to  the  true 
spectrum. 

An  electronic  version  of  this  work  in  its  full  depth  is  included  here.  The  same  paper 
can  be  found  in  the  appendix  of  the  printed  report  and  is  entitled:  Spectral  balance:  A 
frequency  domain  framework  for  analysis  of  nonlinear  dynamical  systems. 
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3.5  Convergence 

Robust  design  of  complex,  interconnected  systems  requires  efficient  computational 
methods  for  uncertainty  propagation.  Knowledge  of  the  probability  distribution  function 
(pdf)  of  key  output  variables  derived  from  the  knowledge  of  input  variables  enables 
better  decision  making  during  design.  Existing  methods  for  propagation  of  uncertainty 
such  as  Monte-Carlo,  polynomial  chaos  and  stochastic  surface  response  methods  scale 
poorly  with  problem  size  and  are  computationally  very  demanding  for  complex  systems 
[!]•  ' 

In  order  to  overcome  the  barrier  of  computational  effort,  a  new  approach  based  on 
decomposition  of  complex  systems  has  been  emerging  as  a  promising  direction.  In  this 
approach,  a  complex  system  is  first  decomposed  using  graph  theoretic  methods  into 
subsystems  connected  in  series,  parallel  or  feedback  [2],  Then,  a  block-by-block 
propagation  is  performed  accounting  for  the  dependencies  among  variables  [3].  Such  a 
method  exploits  the  underlying  hierarchical  structure  of  a  complex  system  and  provides 
the  flexibility  to  use  different  methods  for  different  subsystems  of  the  original  system. 

Feedback  loops  at  the  system  level,  or  encompassing  a  large  number  of  subsystems,  pose 
a  significant  risk  to  the  block-by-block  framework  since  the  structural  decomposition 
methods  described  in  [2]  identify  the  feedback  loops  as  a  single  subsystem.  The 
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decomposition  method  based  on  the  Jacobian  value  and  the  associated  graph  Laplacian, 
also  described  in  [2],  can  identify  weak  feedback  connections  only.  However,  moderate 
to  strong  feedback  due  to  a  controller  is  common  in  many  engineering  systems.  The 
ability  to  propagate  uncertainty  through  such  feedback  loops  block-by-block  without 
having  to  solve  the  closed  loop  system  in  entirety  is  very  attractive.  Thus,  there  is  a  need 
for  alternative  approaches  to  speeding  up  the  computation  for  uncertainty  propagation  in 
feedback  systems. 

In  this  work,  we  developed  an  iterative  method  to  compute  the  probability  density 
function  of  the  output  of  a  feedback  loop  from  the  probability  density  function  of  the 
input.  We  derived  the  problem  setup  and  formulated  the  iteration  equations  by  abstracting 
a  computational  scheme  conceived  and  implemented  by  Huzmezan  and  Kalmar-Nagy 
[3].  We  proved  the  point- wise  convergence  of  the  iteration  to  the  true  closed  loop 
probability  density  function  under  the  assumption  that  the  loop  operator  was  contractive. 

Additionally,  we  considered  the  case  where  there  is  additional  parametric  uncertainty  in 
the  loop  operator.  We  showed  that  the  problem  could  be  cast  as  a  system  of  iterated 
random  functions.  Based  on  the  results  in  [4,5],  we  claim  that  the  iteration  converges  to 
the  solution  of  the  closed  loop  probability  density  function  under  the  assumption  that  the 
loop  operator  is  contractive  on  an  average.  The  proof  of  this  claim  and  extension  to 
dynamic  systems  will  be  considered  in  future  work. 

An  electronic  version  of  this  work  in  its  full  depth  is  included  here.  The  same  paper 
can  be  found  in  the  appendix  of  the  printed  report  and  is  entitled:  An  Iterative  Method  for 
Propagation  of  Probability  Distributions  in  Feedback  Systems. 
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3.6  Conservative  Systems 


This  work  studied  the  effect  of  uncertainty,  using  random  perturbations,  on  area 
preserving  maps  of  Rn  to  itself.  The  focus  was  on  the  standard  map  and  a  discrete 
Duffing  oscillator  in  R2  as  specific  examples.  We  related  the  level  of  uncertainty  to  the 
large-scale  features  in  the  dynamics  in  a  precise  way.  We  also  studied  the  effect  of  such 
perturbations  on  bifurcations  in  such  maps.  The  main  tools  used  for  these  investigations 
was  a  study  of  the  eigenfunction  and  eigenvalue  structure  of  the  associated  Perron- 
Frobenius  operator  along  with  set  oriented  methods  for  the  numerical  computations. 

Uncertainty  is  obviously  of  central  importance  for  dynamic  and  control  systems 
both  from  a  theoretical  and  a  practical  point  of  view.  While  this  topic  has  received  much 
attention  in  the  control  community,  there  is  still  much  to  learn  about  the  effect  of  random 
perturbations  on  such  systems.  This  work  focused  our  attention  on  perturbations  of 
conservative  system  as  a  first  step  towards  studying  mechanical  systems  (including 
molecular  systems)  in  the  presence  of  uncertainty.  An  important  aspect  of  our  approach 
was  to  try  to  extract  the  key  dynamical  features  that  survive  in  a  noisy  environment.  One 
can  make  the  case  that  such  features  are  the  most  important  ones  to  compute. 

An  electronic  version  of  this  work  in  its  full  depth  is  included  here.  The  same  paper 
can  be  found  in  the  appendix  of  the  printed  report  and  is  entitled:  Uncertainty  in  the 
Dynamics  of  Conservative  Maps. 


This  work  showed  that  one  could  effectively  use  the  theory  and  computation 
associated  with  the  Perron-Frobenius  operator  to  study  the  effect  of  uncertainty  on  area 
preserving  maps  and  illustrated  the  methods  on  the  standard  map  and  the  discrete  Duffing 
oscillator.  The  level  of  uncertainty  was  related  in  a  quantitative  way  to  the  large-scale 
features,  often  ones  that  are  the  most  important  to  compute.  This  work  also  studied  the 
way  in  which  uncertainty  affects  the  bifurcation  of  such  maps. 


3.7  Applications 


3.7.1  DNA 

We  investigated  the  dynamics  of  DNA  molecules  using  a  coarse-grained  model 
presented  in  [1].  In  particular  we  analyzed  the  robustness  of  the  beneficial  DNA 
dynamics  with  respect  to  uncertainty  in  the  model  parameters  and  the  initial  conditions.  It 
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was  found  that  the  DNA  model  shows  a  remarkable  robustness  of  the  beneficial 
dynamics  to  uncertainty  in  model  parameters  and  the  initial  conditions.  The  beneficial 
behavior  enables  coding  of  the  genetic  material  and  involves  a  synchronous  behavior  of 
the  molecules. 

Reference  for  this  section: 

[1]  Yakushevich,  L.V.,  Savin,  A.V.,  and  Manevitch,  L.I.,  “Nonlinear  dynamics  of  topological  solitons 
in  DNA,”  Physical  Review  E  66,  016614 ,  2002. 

3.7.2  Molecular  Models 

We  summarize  the  efforts  made  in  the  following  three  areas  of  molecular  modeling 
in  the  context  of  uncertainty  analysis: 

(a)  We  have  identified  systems  in  molecular  mechanics  for  conformational  analysis  by 
classical  molecular  dynamics.  The  globally  optimum  configuration  is  subject  to 
uncertainty  in  the  inter-atomic  potential.  We  propose  to  test  polynomial  chaos  on 
classical  molecular  dynamics  calculations  on  the  nitrogen  molecule  and  the  more 
complicated  trimer  molecule  CS2. 

(b)  We  have  been  studying  the  quantification  of  uncertainty  in  classical  and  quantum 
molecular  dynamics  theories  in  relation  to  property  predictions  in  molecular  systems. 
These  theories  involve  the  simulation  of  dynamical  systems  in  time  that  may  cover 
such  processes  as  bond-breaking,  bond-formation,  chemical  reaction,  and  diffusion. 
These  theories  are  general  but  our  ultimate  application  targets  are  in  such  areas  as 
proton  diffusion  in  PEM  fuel  cell  membranes  and  hydrogen  absorption  in  metal- 
hydride  materials  for  hydrogen  storage.  A  fundamental  issue  is  uncertainty 
propagation  in  these  dynamical  systems.  The  uncertainty  arises  in  parameters 
pertaining  to  the  semi-empirical  force-fields  if  classical  molecular  dynamics  are  used; 
this  can  affect  the  accuracy  of  transport  properties,  such  as  diffusivity,  estimated  by 
molecular  dynamics.  We  have  identified  this  issue  as  a  first  problem  to  be  rigorously 
analyzed  by  polynomial  chaos  methods.  We  plan  to  compute  the  errors  in  the 
diffusivity  obtained  from  classical  molecular  dynamics  calculations  subject  to  random 
errors  in  the  parameters  comprising  the  functional  form  of  the  inter-atomic  potential. 
We  plan  to  start  with  a  simple  Lennard- Jones  system  with  uncertainty  in  the  collision 
diameters  and  the  well-  depths  and  propose  to  extend  the  study  to  cover  more 
sophisticated  potentials  for  metal-ceramics  of  interest  to  UTC  and  solid-oxide  fuel 
cells  within  the  context  of  polynomial  chaos. 

(c)  When  molecular  modeling  calculations  are  done  at  a  more  fundamental  level  (e.g., 
Hohenberg-Kohn  density  functional  theory  at  the  electronic  level)  there  is  the 
uncertainty  in  the  form  of  the  exchange-  correlation  functional.  Despite  forty  years  of 
density  functional  theory,  the  local  density  approximation  and  the  generalized 
gradient  approximations  are  the  only  parameterizations  available  for  the  exchange 
correlation  functional.  We  know  that  these  are  not  satisfactory  for  instance  in  the 
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accurate  calculations  of  cohesive  energies,  reaction  barrier  heights  etc.  The 
application  of  polynomial  chaos  to  study  uncertainty  propagation  due  to  the  exchange 
correlation  functional  is  an  excellent  idea.  But,  it  cannot  be  done  without  an 
understanding  of  the  uncertainty  in  the  density  functional.  We  propose  to  consider 
using  the  semi-empirical  density  functionals  that  the  chemists  have  constructed  by 
fitting  a  guessed  at  form  with  parameters  to  the  results  of  configuration-interaction 
computations  for  small  molecular  systems.  This  can  be  taken  as  a  standard  for 
comparisons.  Then  the  various  physically  motivated  functionals  would  deviate  from 
that  standard,  and  one  could  study  the  propagation  of  those  deviations  through  the 
Car-Parrinello  quantum  dynamics  computation.  This  area  is  of  interest  to  UTC  as 
Car-Parrinello  calculations  are  used  in  hydrogen  storage  calculations  and  in  the 
calculation  of  chromia  ion  migration  in  the  cathodes  of  solid  oxide  fuel  cells. 

An  electronic  version  of  this  work  in  its  full  depth  is  included  here.  The  same  paper 
can  be  found  in  the  appendix  of  the  printed  report  and  is  entitled:  Dynamics  in  Molecular 
Modelling  and  the  Scope  for  Uncertainty  Analysis  by  Dellnitz-Preis  and  Polynomial 
Chaos  Methods. 


3.7.3  Aeroacoustics 

For  this  application  we  studied  a  thermo-acoustic  model  on  a  cylindrical  or  annular 
geometry,  capable  of  modeling  instabilities  of  tangential  acoustic  modes.  The  model 
accounts  for  nonuniform  density,  damping,  rotational  flow,  and  heat-release  coupling.  It 
is  shown  that  deliberately  introducing  spatial  variations  in  some  quantities  has  a  similar 
effect  to  adding  damping  to  the  system.  The  effects  of  these  symmetry-breaking  concepts 
are  evaluated  on  the  model  through  linear  analysis  and  the  net  amount  of  additional 
damping  is  computed.  We  show  how  various  symmetry-breaking  concepts  are  robust 
with  respect  to  the  uncertainty  in  the  model  parameters  and  we  examine  propagation  of 
uncertainty  with  respect  to  a  recently  defined  measure  of  uncertainty. 

Thermo-acoustic  instabilities  in  gas  turbine  and  rocket  engines  develop  when 
acoustic  waves  in  combustors  couple  with  an  unsteady  heat-release  field  in  a  positive 
feedback  loop.  For  a  summary  of  active  control  of  thermo-acoustic  instabilities  see  [3]. 
Thermo-acoustic  modeling  and  control  is  well-studied  for  axially  extended  combustion 
chambers,  as  in  [5],  [7],  [14],  [10],  where  the  acoustic  to  heat  release  coupling  is 
dominated  by  longitudinal  acoustic  modes.  However,  comparably  less  attention  has 
focused  on  thermo-acoustic  modeling  in  combustion  chambers  with  annular  or  cylindrical 
geometries.  Our  efforts  for  the  aeroacoustic  application  focused  on  the  development  a 
low-order  thermo-acoustic  model  on  a  circular  geometry,  similar  to  that  of  [1]. 
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Often  thermo-acoustic  instabilities  are  dominated  by  a  few  natural  acoustic  modes, 
which  can  accurately  be  modeled  with  a  low-dimensional  model.  Accurate  heat-release 
models  are  difficult  and  time-consuming  to  implement.  A  low-order  thermo-acoustic 
model,  properly  calibrated  with  acoustic  data,  can  provide  insight  into  the  possibly 
deleterious  acoustic-heat-release  coupling  and  may  provide  a  platform  for  fast  evaluation 
of  preliminary  design  concepts. 

We  began  be  defining  symmetry-breaking  as  the  deliberate  introduction  of  spatial 
variations  in  the  system  parameters  in  order  to  change  the  stability  properties.  Recent 
work  has  focused  on  analysis  of  heterogeneous  distributed  systems  [6],[9],[11]. 
Symmetry-breaking  is  commonly  referred  to  as  mistuning  in  the  literature  regarding  the 
dynamics  of  arrays  of  turbine  blades  on  a  disk.  Studies  of  stability  properties  of  turbine 
blade  flutter  through  the  introduction  of  spatial  nonuniformities  has  appeared  in  [2],  [16]. 
Optimal  mistuning  in  arrays  of  bladed  disks  has  appeared  in  [15],  [17].  A  study  of  the 
effects  of  asymmetry  on  compressor  stall  inception  has  appeared  in  [8]. 

As  in  the  case  of  mistuning  in  arrays  of  bladed  disks  in  turbines,  this  form  of 
passive  control  is  often  more  feasible  than  implementing  an  active  control  scheme.  This 
may  also  be  true  for  the  case  in  combustion  chambers,  where  high  temperatures  prohibit 
adequate  sensing  and  may  damage  the  actuators  required  for  active  control.  Furthermore, 
We  have  found  that  symmetry-breaking  can  be  a  more  cost-effective  means  of  stability 
enhancement. 

An  electronic  version  of  this  work  in  its  full  depth  is  included  here.  The  same  paper 
can  be  found  in  the  appendix  of  the  printed  report  and  is  entitled:  Symmetry-Breaking 
and  Uncertainty  Propagation  in  a  Reduced  Order  Thermo-acoustic  Model. 
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(c)  Transitions 

The  symmetry  breaking  concepts  to  reduce  sensitivity  of 
dynamical  systems  to  uncertainty  developed  under  this  contract  were 
applied  in  two  internally  funded  projects  to  investigate  design  of 
aeroengine  combustors  robust  to  thermoacoustic  instabilities. 

Using  advanced  mathematical  concepts  developed  for  uncertainty 
propagation  through  dynamical  systems  UTRC  and  PW  are  currently 
looking  at  validation  and  verification  (V&V)  of  model  based  controllers 
for  life  extension  of  hybrid  systems  such  as  gas  turbines. 
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Abstract — Physics  based  models  are  often  converted 
to  monolithic  systems  of  uncertain  nonlinear  differen- 
tial/algebraic  equations.  Graph  decomposition  methods  can 
be  used  to  decompose  such  system  into  subsystems  evolv¬ 
ing  on  different  time  scales.  This  time  scale  separation 
can  be  exploited  to  increase  computational  efficiency  when 
propagating  input  uncertainty  in  a  subsystem-by-subsystem 
manner.  In  this  paper  the  propagation  of  uncertain  inputs 
through  series,  parallel  and  feedback  interconnections  of  dy¬ 
namical  systems  with  simple  asymptotic  behavior  is  studied 
by  employing  discrete  density  mapping  (analogous  to  the 
input-output  Perron-Frobenius  operator).  A  process  con¬ 
trol  example  is  used  to  illustrate  the  method. 

Index  Terms — Perron-Frobenius,  Density  Mapping,  Un¬ 
certainty  Propagation,  Networks  of  Dynamical  Components 


I.  Introduction 

Enhancing  engineering  performance  and  productivity 
through  systematic  and  integrated  designs  that  account 
for  the  effects  of  uncertainty  is  a  key  economic  prior¬ 
ity  [1].  A  unified  paradigm,  the  ’’analytical  systems  en¬ 
gineering”,  recently  formulated  at  the  United  Technolo¬ 
gies  Research  Center,  has  identified  uncertainty  propaga¬ 
tion  through  networks  of  nonlinear  components  as  an  es¬ 
sential  component.  This  view  is  also  echoed  in  [2]  or  in 
[3],  where  large-scale,  interconnected  systems  arc  designed 
using  model-based  techniques  that  employ  uncertainty  de¬ 
scriptions  explicitly.  Since  all  system  models  have  varying 
levels  of  uncertainty  [4],  designers  often  use  large  safety 
margins,  which  result  in  more  complex  and  expensive  sys¬ 
tems  [5].  As  a  result,  a  natural  path  in  modern  systems 
design  is  to  make  decisions  on  the  best  system  structure 
from  the  perspective  of  greater  robustness  to  uncertainty. 

Before  moving  towards  design  addressing  uncertainty 
analysis  in  a  model-based  complex  engineering  systems 
framework  is  important.  This  has  been  traditionally  been 
addressed  by  Monte  Carlo  (MC)  like  methods.  This  classic 
approach  (see  [6])  employs  a  large  number  of  simulations 
with  a  random  selection  of  variables  from  their  prescribed 
distribution  (parametric  or  empirical). 

Unfortunately,  uncertainty  propagation  techniques  using 
MC  methods,  even  in  their  advanced  forms,  do  not  scale 
well  with  system  size.  The  natural  choice  in  such  cases  is 
to  break  the  large  system  into  pieces.  More  precisely,  the 
physics  based  models,  which  are  often  converted  to  mono- 
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lithic  systems  of  uncertain  nonlinear  differential/algebraic 
equations,  are  to  be  decomposed  using  graph  decomposi¬ 
tion  methods  into  subsystems  evolving  on  different  time 
scales,  as  mentioned  in  [7]. 

In  the  authors’  view  the  time  scale  separation  can  be 
exploited  to  increase  computational  efficiency  when  prop¬ 
agating  input  uncertainty  in  a  subsystem-by-subsystem 
manner.  This  approach  has  been  also  advocated  in  [8], 
where  arbitrary  interconnections  of  multivariable  systems 
(represented  either  in  a  continuous  or  discrete  form)  with 
nonlinear  or  linear  dynamics  (nonlinear  time  varying,  dis¬ 
tributed  linear  time  invariant  or  lumped  linear  time  invari¬ 
ant)  are  decomposed  into  aggregate,  strongly  connected 
subsystems.  Following  this  procedure  each  subsystem  is 
addressed  using  ’’the  minimum  set  technique”  and  trans¬ 
formed  into  a  typical  feedback  interconnection. 

Alternatively,  in  [9],  a  Recursive  Projection  Method 
(RPM)  is  developed  to  solve  nonlinear  parameter  problems 
for  which  the  convergence  is  achieved  for  certain  parame¬ 
ter  values.  The  correction  applied  for  divergent  domains 
involves  Newton’s  integration  method.  RPM  provides  reli¬ 
able  results,  when  the  number  of  divergent  modes  is  small 
comparative  to  the  system’s  dimension.  The  algorithm  has 
been  successfully  demonstrated  on  folds,  bifurcations  and 
unstable  system  branches.  This  approach  is  believed  by  its 
authors  to  greatly  accelerate  iteration  convergence. 

Along  the  same  thinking  path,  extensions  from  the  back¬ 
ward  Euler  formula,  conventionally  used  to  obtain  a  system 
of  nonlinear  algebraic  equations  from  an  original  system  of 
nonlinear  algebraic  differential  equations,  can  be  grouped 
under  the  waveform  relaxation  (WR)  method.  This  tech¬ 
nique  has  been  successfully  employed  in  [10]  to  address 
large  scale  systems,  such  as  integrated  circuits.  The  iter¬ 
ative  WR  method  decomposes  the  system  in  several  dy¬ 
namical  subsystems,  which  arc  independently  analyzed  for 
the  entire  time  interval.  This  method  comes  with  sufficient 
convergence  guarantees,  also  revealed  in  [10]. 

In  [11]  graph  theory  is  also  used  in  the  context  of  au- 
tocatalytic  networks/sets  to  classify  the  uncertainty  of  the 
network  and  predict  its  influence  over  short  and  medium 
time-scales.  Therefore,  looking  at  how  networks  evolve 
with  time,  more  precisely  looking  at  their  dynamics,  it  can 
be  beneficial  from  a  computation  speed  perspective.  It  is 
essential,  when  dealing  with  irreducible  (i.e.  strongly  con¬ 
nected)  graphs  (e.g.  dynamical  subsystems),  to  perform 
related  computations  in  corresponding  time  scales  [12]. 

Finding  a  subset  of  the  state  space  of  a  dynamical  system 
where  typical  trajectories  stay  longer  before  entering  differ¬ 
ent  regions,  conventionally  called  almost  invariant  sets  [13], 
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entitles  a  macroscopic  behavior  analysis  of  dynamical  sys¬ 
tems  networks.  In  effect  this  separation  permits  individual 
sub-system  uncertainty  propagation  studies.  Algorithms 
for  containment  of  such  almost  invariant  sets  are  presented 
in  [13]. 

To  address  slow  convergence  rate  and  clustering  issues 
associated  with  Monte  Carlo  methods,  new  methods  such 
as:  polynomial  chaos  [14]  (i.e.  stochastic  finite  elements), 
stochastic  surface  response  methods  [15]  and  probabilistic 
collocation  methods  [16]  were  used  with  significant  suc¬ 
cess.  Complementing  these  approaches,  the  propagation 
of  uncertainty  in  the  distribution  of  initial  conditions  for  a 
network  of  dynamical  systems  can  be  studied  in  an  exact 
manner  using  the  Liouville’s  equation  as  in  [17]. 

A  bottle  neck  in  the  application  of  uncertainty  propaga¬ 
tion  methods  in  a  subsystem-by-subsystem  way  is  associ¬ 
ated  with  the  existence  of  correlated  signals  due  to  typical 
parallel  or  feedback  connections  found  in  conventional  sub¬ 
systems.  To  address  this  issue,  in  this  paper,  the  propa¬ 
gation  of  uncertain  inputs  through  interconnections  of  dy¬ 
namical  systems,  with  simple  asymptotic  behavior,  is  stud¬ 
ied.  The  method  of  choice  is  a  discrete  density  mapping, 
analogous  to  the  input-output  Perron-Frobenius  operator. 

The  advocated  method  assumes  that  a  large  system 
is  broken  into  components,  which  evolve  on  different 
timescales  or  have  simple  attractor  structure.  Then,  uncer¬ 
tainty  propagation  methods  are  used  to  map  input  distri¬ 
bution  through  components  and  arrive  to  a  stationary  den¬ 
sity  for  the  states/outputs  of  interest.  In  a  chain  topology 
only  the  need  for  a  transfer  operator  is  observed.  This  is  in 
contrast  with  the  parallel  connection,  which  adds  complex¬ 
ity  through  the  necessity  of  summing  correlated  signals.  In 
the  case  of  feedback  structures  this  complexity  is  further 
augmented  by  the  requirement  for  converge  of  the  itera¬ 
tions. 

Following  the  introduction,  Section  II  makes  reference 
to  chains,  parallel  and  feedback  connections  which  arc  ap¬ 
proached  with  the  density  mapping  method  described  in 
Section  III.  In  Section  III  B  convergence  guarantees  are 
addressed.  To  illustrate  the  method  an  industrial  process 
control  example  abstraction  is  offered  in  Section  IV.  The 
conclusions  are  stated  in  Section  V. 

II.  Signal  Flow  Decomposition:  Chains,  Parallel 
and  Feedback  Connections 

The  starting  assumption  for  the  uncertainty  propaga¬ 
tion  method  application  is  that  a  nonlinear  system  under 
investigation  is  broken  into  interconnected  components. 

In  this  case  the  simplest  topology  that  can  arise  is  a 
chain,  as  shown  in  Figure  1.  In  this  case  the  output  of  a 
block  serves  as  input  to  the  following  block.  If  the  input- 
output  maps  fi,  ...  are  known  for  all  blocks,  then 
the  output  of  the  chain  can  simply  be  calculated  as  the 
composition  of  these  maps  f\of2o...  acting  on  the  input  of 
the  first  block.  To  consider  uncertain  inputs,  it  is  necessary 
to  extend  the  notion  of  single  input-output  mapping  to 
probability  densities.  The  resulting  formalism,  analogous 
to  the  Perron-Frobenius  operator  is  discussed  below. 


Fig.  1.  Density  propagation  through  chains 
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Fig.  2.  Parallel  summation 


The  computational  advantage  of  propagating  input  den¬ 
sities  in  a  block-by-block  manner,  which  corresponds  to  the 
composition  of  maps,  becomes  clear  when  the  dynamics 
of  different  blocks  include  completely  different  timescales. 
Complications  arise  when  considering  a  parallel  connec¬ 
tion,  see  Figure  2.  In  this  case,  the  resulting  probability 
densities  (output  uncertainties)  have  to  be  summed.  For 
densities  of  independent  random  variables,  this  operation 
would  be  a  simple  convolution.  However,  when  these  densi¬ 
ties  are  dependent,  correlation  information  should  be  used 
to  sum  them  correctly,  fact  also  discussed  in  the  next  sec¬ 
tion. 


Fig.  3.  Feedback  loop  with  loop  operator  G 
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Finally,  the  feedback  connection,  shown  in  Figure  3,  can 
be  thought  of  as  infinite  series  of  parallel  connections  for 
which  the  convergence  aspects  are  discussed  in  Subsec¬ 
tion  B. 


III.  The  Density  Mapping  Method 

The  paper  proposes  a  simple  algorithm  that  can  be  used 
to  propagate  input  uncertainty  through  nonlinear  compo- 
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nents.  While  this  algorithm  can  be  thought  of  as  a  simple 
implementation  of  the  Perron-Frobenius  operator  [18],  it 
has  been  extended  to  account  for  the  mapping  of  the  dis¬ 
tribution  support  (akin  to  ccll-to-cell  mapping  [19]).  This 
extension  was  necessary  to  enable  operations  with  corre¬ 
lated  densities  (e.g.  parallel  or  feedback  loop  connections). 

To  understand  the  proposed  technique  first  we  review 
the  Perron- Frobcnius  operator  of  a  scalar  map  in  the  con¬ 
text  of  the  proposed  input-output  uncertainty  mapping. 
The  one-dimensional  map: 

Un+ 1  =  /  («n)  (1) 

corresponds  to  the  action  of  a  block  upon  un ,  the  system 
input.  This  action  results  by  applying  on  un  the  mapping 
corresponding  to  the  asymptotic  dynamics  of  the  system  / 
and  has  as  result  un+ i,  the  output.  The  uncertainty  of  un 
will  be  represented  by  its  probability  density  G.  Therefore 
the  map  has  a  corresponding  Perron-Frobenius  operator  U 
acting  on  G : 

Gn+1  =UGn  —  j 8  (u  —  f  (x))  Gn  (x)  dx  = 

xp  Gn  (fa  1  (u))  /<y\ 

u 

where  the  summation  should  be  taken  over  all  inverse 
branches  f~l  ( u ). 


Fig.  4.  Density  mapping 

In  the  following  a  straightforward  numerical  implemen¬ 
tation  of  this  mapping  is  described.  Taking  a  cell-based 
approximation  of  the  input  probability  density  (histogram) 
{(r;,7v+i) ,  Gj}  with  cell-width  w  a  collection  of  rectangles 
(bins)  is  produced.  Consecutively  the  mapping  /  is  ap¬ 
plied  to  the  corner  points  of  all  rectangles.  Their  resulting 
heights  are  derived  from  the  density  conservation  condi¬ 
tion.  This  procedure  yields  a  new  set 

«/(n),/(n+.)).M  '■»=|/(T,)a;(r,tl)|  W 

If  the  map  /  is  many-to-one,  this  collection  will  contain 
overlapping  rectangles. 

The  overall  output  density  can  be  produced  by  ’re- 
binning5,  a  procedure  which  involves  finding  the  total  area 


over  a  bin  on  the  support  of  the  resulting  distribution 
[min  /  (ri) ,  max  /  (t*)].  This  corresponds  to  summing  over 
the  inverse  branches  of  the  map.  Re-binning  is  useful  when 
propagating  densities  in  a  system  with  chain  topology.  On 
the  other  hand,  re-binning  destroys  information  about  the 
mapping  of  the  original  support  into  the  support  of  the  out¬ 
put  in  the  same  fashion  as  polynomial  chaos  methods.  This 
information  is  crucial  when  dealing  with  densities  mapped 
through  parallel  or  feedback  connections. 

A.  Summing  Correlated  Probability  Densities 

Consider  the  parallel  connection  shown  in  Fig¬ 
ure  2.  At  the  summing  junction  we  have  two  lists 
{(fi(xi),fi{xi+1)),hi},  {(f2{xi),f2{xi+l)),hi}  pro- 

duced  from  the  same  input  density  G.  Note  that  the  lists 
contain  the  same  /ij’s,  because  of  the  density  conserva¬ 
tion.  Consequently,  to  produce  the  parallel  structure  out¬ 
put  density,  the  sum  of  these  lists  are  then  taken  as 

{(/i  (z*) ,  fi  (zi+i))  ,hi}  ©  {(/2  (Zi) ,  /2  (Zi+i))  ,hi}  = 

=  {(/i  (zi)  +  h  (Zi) ,  h  (z,+i)  +  h  (zi+i)) ,  hi}  (4) 

This  represents  the  list  produced  by  the  simple  application 
of  the  operator  fi  -f  /2  to  the  original  density. 

B.  Feedback  Loop:  Convergence  Issues 

As  discussed  the  feedback  connection  represents  a  natu¬ 
ral  extension  of  dealing  with  the  summation  of  correlated 
signal  density  functions.  Obtaining  the  feedback  signals 
involves  a  summation  of  a  set-point,  generally,  with  a  pro¬ 
cessed  feedback  signal  that  has  been  generated  as  a  result 
of  that  set  point.  Propagating  densities  instead  of  signals 
through  the  loop  is  facilitated  by  keeping  track  of  the  un¬ 
certainty  propagation  in  terms  of  density  lists.  For  com¬ 
putation  purposes  a  finite  number  of  iterations  is  expected 
before  convergence  is  realized.  To  minimize  the  compu¬ 
tational  load  management  of  the  bin  size  for  the  original 
distribution  is  required.  Under  such  conditions  the  sav¬ 
ings  obtained  arc  significant,  i.e.  two  orders  of  magnitude 
smaller  versus  traditional  methods  such  as  Monte  Carlo. 

A  rigorous  proof  for  the  convergence  of  the  iterative 
scheme  proposed  is  offered  in  [20].  For  completeness  a 
sketch  of  the  proof  is  presented  here  for  feedback  systems 
falling  under  the  class  presented  in  Figure  3  for  which  u  is 
the  input,  x  the  internal  signal  and  y  =  fx,  where  /  is  an 
operator,  the  output. 

The  above  signals  have  associated  corresponding  prob¬ 
ability  density  functions  (pdf)  (i.e.  for  u  the  pdf  is  Gu). 
Applying  the  Perron-Frobenius  theorem  on  the  closed  loop 
relation  x  —  (I  -  f)~xu,  the  pdf  of  z,  Gx,  can  be  written 
as: 

Gx  =  GU((I  -  f)x)\I  —  df\  (5) 

where  df  is  the  derivative  of  /  and  |  •  |  denotes  the  absolute 
value  of  the  determinant.  The  proof  developed  in  [20] 
shows  that  a  sequence  iteration  for  Gx  converges  to  the 
closed  loop  expression  in  Eq  5,  which  is  the  correct  closed 
loop  solution. 
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Fig.  5.  Typical  process  control  structure  and  its  corresponding  asymptotic  map 


Note  that  the  operator  /  can  be  nonlinear.  Also  note 
that,  as  presented  in  [20],  the  multidimensional  generaliza¬ 
tion  of  this  computational  method  is  possible. 

C.  Benefits  of  a  Block-by-block  Propagation  in  a  Multi 
Time  Scale  System 

Another  benefit  of  using  the  proposed  approach  is  re¬ 
vealed  in  a  multi  time  scale  system.  For  such  n-equation 
system  a  simple  calculation  can  quantify  the  computational 
benefits. 

Assuming  that  to  compute  the  solution  of  one  equation 
requires  (in  average)  k  operations  and  assuming  that  the 
step-size  required  to  resolve  the  fastest  timescale  is  dtfast 
it  results  that  to  integrate  the  system  to  tfinai  requires 
approximately  n  x  k  x  steps.  Further,  assuming 

an  overall  system  decomposition  into  10  subsystems,  with 
fastest  timescale  dtjast  separated.  The  number  of  oper¬ 
ations  required  to  integrate  this  system  is  approximately 
■—xnxkx  hfy—t ,  which  represents  an  order  of  magnitude 
computation  speed  up. 

IV.  Process  Control  Example 

Well  acclaimed  techniques  steaming  from  robust  control 
theory  such  as  loop-shaping  or  p  analysis  and  synthe¬ 
sis  arc  considering,  in  general,  worst  case  scenarios  and  do 
not  use,  most  of  the  time,  existent  additional  information 
(i.e.  probabilistic  knowledge)  about  uncertainty.  As  shown 
in  this  paper,  a  direct  characterization  of  uncertainty  is  a 
possible  and  useful  alternative.  Using  this  approach,  for  a 
control  systems  framework,  is  a  useful  exercise,  which  can 
be  further  generalized  to  complex  networks  of  dynamical 
components.  The  example  for  uncertainty  propagation  has 
been  generated  through  the  abstraction  of  a  typical  indus¬ 


trial  controls  closed  loop.  This  process  control  structure, 
see  Figure  5,  involves:  inner  control  loops,  feedforwards 
and  outer  control  loops. 

Using  the  developed  tools  to  understand  how  input 
uncertainty  propagates  through  this  network  of  dynami¬ 
cal  components  requires  each  component’s  nonlinear  map. 
These  maps  can  be  extracted  from  the  asymptotic  behavior 
of  the  individual  nonlinear  dynamic  subsystem. 

Therefore,  the  network  structure,  also  displayed  in  Fig¬ 
ure  5,  is  captured  through  steady  state  models  displaying 
series,  parallel  and  feedback  interconnections  of  nonlinear 
components. 

Recall  that  the  density  mapping  method  operates  in 
terms  of  density  lists,  which  once  mapped  through  a  block 
are  constructing  other  lists.  This  procedure  replaces  tradi¬ 
tional  signal  propagation  through  blocks.  The  interconnec¬ 
tions  of  the  overall  system,  first  decomposed  in  standard 
building  blocks,  arc  represented  with  in  a  similar  fashion 
as  in  the  Matlab  Control  Systems  Toolbox. 

The  evolution  of  an  initial  uncertain  input  distribution 
(solid  thick  line),  presented  in  Figure  6,  localized  at  outer 
loop  set  point  level  is  mapped  in  a  successive  manner 
through  the  network  components.  The  result  showing  the 
successive  convergence  to  the  known  closed  loop  solution 
(also  solid  thick  line),  known  in  this  simulation  example,  is 
displayed  in  Figure  6  (dotted  lines).  Note  that  to  produce 
this  plot  the  density  lists  were  ’re-binned’  as  presented  in 
Section  III,  based  on  the  same  support  and  bin  size  as  for 
the  original  distribution. 

V.  Conclusions 

The  propagation  of  uncertain  inputs  through  scries,  par¬ 
allel  and  feedback  interconnections  of  dynamical  systems 
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Fig.  6.  Density  mapping  iterations  convergence 

with  simple  asymptotic  behavior  has  been  studied  by  em¬ 
ploying  the  discrete  density  mapping  method.  Conver¬ 
gence  guarantees  further  explored  in  [20]  reflect  the  prac¬ 
tical  observations  made  when  developing  the  code. 

The  proposed  method  achieves  at  least  an  order  of  mag¬ 
nitude  computational  gains  versus  Monte  Carlo  methods. 
Understanding  the  role  of  density  lists  in  dealing  with  cor¬ 
related  signal  density  functions  provides  a  framework  for 
an  extension  of  current  work  to  integrate  PC  methods. 
The  basic  idea  of  this  future  approach  is  to  decompose  the 
distribution  of  an  uncertain  input  into  small  pieces  of  uni¬ 
form  distributions  and  propagate  these  distributions  based 
on  the  PC  methods.  In  this  way  the  correlation  informa¬ 
tion  is  preserved  for  later  operation  on  dependent  densi¬ 
ties.  Given  the  uniform  distribution  characteristic  of  every 
cell  the  PC  expansion  will  be  performed  on  an  orthonor¬ 
mal  Legendre  basis.  This  extension  opens  the  road  for  the 
propagation  of  uncertain  inputs  through  interconnections 
of  uncertain  dynamical  systems  with  uncertain  initial  con¬ 
ditions. 

Extensions  of  the  current  method  to  multivariable  non¬ 
linear  dynamical  systems  are  possible.  Considering  a  sin¬ 
gle  multivariable  nonlinear  block  note  that  the  uncertainty 
propagation  method  presented  generalizes  to  volumes  in¬ 
stead  of  areas  since  the  underlaying  principle  is  the  con¬ 
servation  of  probability.  Also  dealing  with  products  is 
achieved  by  using  a  logarithmic  conversion  followed  by  the 
application  of  knowledge  acquired  in  the  case  of  summing 
correlated  signal  density  functions. 
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Abstract — Uncertainty  propagation  in  complex,  interconnected 
dynamical  systems  systems  can  be  performed  more  efficiently 
by  decomposing  the  network  based  on  the  hierarchy  and/or 
the  strength  of  coupling.  In  this  paper,  we  first  present  a 
structural  decomposition  method  that  identifies  the  hierarchy  of 
subsystems.  We  briefly  review  the  notion  of  horizontal-vertical 
decomposition  (HVD)  or  strongly  connected  components  (SCC) 
decomposition  of  a  dynamical  system  and  describe  algorithms 
based  on  Markov  chain  theory  and  graph  theory  to  obtain 
the  HVD  from  the  equation  graph  of  the  system.  We  also 
present  a  non-structural  decomposition  method  to  identify  the 
weakly  connected  subsystems  of  a  system  based  on  the  Laplacian 
of  a  graph  derived  from  the  Jacobian.  While  most  of  prior 
efforts  in  this  direction  concentrated  on  stability,  robustness  and 
concrete  results  were  limited  to  linear  systems,  we  use  it  for 
uncertainty  propagation  and  study  of  asymptotic  behavior  of 
nonlinear  interconnected  systems.  We  illustrate  the  two  methods 
using  a  fuel  cell  system  example.  These  two  methods  provide  a 
framework  for  efficient  propagation  of  uncertainty  in  complex 
nonlinear  systems. 

I.  Introduction 

Complex  physical  systems  often  possess  an  inherent  struc¬ 
ture  from  which  they  derive  the  robustness  of  their  behavior. 
For  example,  the  DNA  molecule  has  a  strong  backbone 
structure  that  allows  global  modes  [1]  and  gene  regulatory 
networks  have  a  hierarchical  structure  that  allows  complex 
gene  expressions  [2].  Similarly,  complex  engineering  systems 
such  as  coordinated  group  of  vehicles,  sensor  and  communi¬ 
cation  networks  are  designed  by  interconnecting  subsystems 
in  a  systematic  manner. 

There  is  significant  interest  in  robust  design  methods  for 
complex,  interconnected  dynamical  systems.  The  design  of 
such  systems  can  be  significantly  accelerated  by  using  models 
for  robustness  assessment  and  redesign.  This  requires  un¬ 
derstanding  how  uncertainty  in  various  parameters  and  the 
system  model  itself  affects  the  model  outputs.  Computational 
methods  such  as  Monte-Carlo,  polynomial  chaos  and  various 
modifications  of  these  can  be  used  to  propagate  probabilistic 
information  from  the  inputs/parameters  to  the  outputs  of  a 
static  or  dynamical  system.  Unfortunately,  the  computational 
effort  for  these  methods  scales  poorly  with  system  complex¬ 
ity  and  there  is  a  strong  need  for  efficient  alternatives.  A 
promising  direction  is  the  decomposition  of  a  complex  system 
using  graph  theoretic  methods  followed  by  block-by-block 
propagation  with  iterations  when  necessary. 

Graph  theoretic  methods  are  used  widely  in  control  theory 
[3]  computer  science  and  network  systems.  More  recently, 


graph  theoretic  techniques  have  been  employed  in  conjunction 
with  set  oriented  numerical  methods  for  computing  the  almost 
invariant  sets  of  dynamical  systems  [4].  There  is  a  close 
connection  between  Markov  chain  theory  and  graph  theory  [5]. 
We  refer  the  reader  to  standard  text  books  and/or  review  papers 
such  as  [6],  [7]  for  basic  notions  from  algebraic  graph  theory 
and  graph  algorithms  such  as  spanning  and  induced  subgraphs, 
cycles,  connectedness  of  undirected  graphs,  strongly  connect¬ 
edness  and  topological  sorting  of  directed  graphs,  adjacency, 
incidence  and  Laplacian  matrices,  the  significance  of  the 
eigenvalues  of  these  matrices,  reducibility  and  irreducibility 
of  the  adjacency  matrix,  depth  first  search  and  breadth  first 
search. 

The  aim  of  this  paper  is  to  introduce  graph  theoretic 
methods  for  decomposing  a  complex  dynamical  system  into 
hierarchically  and/or  weakly  connected  subsystems  for  the 
purpose  of  uncertainty  propagation  in  systems.  The  structure 
resulting  from  graph  decomposition  with  series,  parallel  or 
feedback  connections  among  subsystems  can  be  exploited  for 
efficient  methods  for  computation  of  uncertainty  propagation 
[8].  We  briefly  review  a  structural  decomposition  algorithm 
[3],  [9]  that  employs  the  so-called  equation  graph  of  a  dynam¬ 
ical  system  in  and  identifies  the  subsystems  and  the  associated 
hierarchy  in  Section  II.  The  notion  of  horizontal-vertical 
decomposition  (HVD)  of  a  dynamical  system  is  introduced 
and  two  algorithms  for  generating  the  HVD  are  presented. 

The  structural  decomposition  method  presented  in  Section  II 
considers  only  the  directionality  of  influence  among  variables 
(/.e.,  whether  a  variable  influences  another  variable)  but  not 
the  strength  of  interaction.  A  limitation  of  this  method,  as 
described  in  Section  II-A,  is  that  it  does  not  yield  any  useful 
structure  for  a  system  or  subsystem  whose  equation  graph  is 
strongly  connected.  However,  this  limitation  can  be  overcome 
using  another  graph  theoretic  method  described  in  Section  III. 
The  latter  method  accounts  for  the  strength  of  coupling  be¬ 
tween  variables  and  identifies  any  weakly  connected  subsys¬ 
tems  using  a  graph  partitioning  technique.  The  two  methods 
can  be  applied  in  conjunction  to  achieve  a  fine  decomposition 
of  a  complex  system.  We  demonstrate  the  methodology  using 
a  reduced  order  fuel  cell  system  model  in  Section  IV. 

II.  Structural  decomposition  methods  for 

DYNAMICAL  SYSTEMS 

Let  (X,fi,T\P)  be  a  dynamical  system  and  J\f  = 
1, ...,  N  C  Z  where  X  -  Xi}  each  X{  a  compact  subset 


of  M.  Let  Hi  :  X  —>  Xi  be  the  usual  projection  to  the  i  —  th 
component. 

In  addition,  \i  =  where  pi  is  a  probabilistic 

measure  on  X{  and  Tl{x,p)  :  X  x  P  — >  X  is  either  a 
discrete-time  family  induced  by  a  map  (in  which  case  t  e  Z) 
or  a  flow  of  a  system  of  first-order  autonomous  system  of 
ordinary  differential  equations.  The  parameters  p  are  defined 
on  P  =  (&ieM  Pi,  where  each  Pi  a  compact  subset  of  E. 
The  change  in  re*  is  given  by  a  function  fi(x),  where  x  €  X. 
This,  of  course  is  the  i  —  th  component  of  either  the  vector 
field  f(x)  —  d/dt\t=o(Tt(x1p))  or  Tl(x,p)  -  x,  in  the  case 
tGZ.  The  Jacobian  matrix  is  defined  as 

J{x,p)  =  Df(x,p), 


or,  in  components, 


Jji(x,p)  = 


where  j  represents  the  column  and  i  the  row  of  the  matrix. 
We  define  the  following  matrix  M(x,p)  from  J{x,p): 


Mji(x,p) 

Mji(x,p) 


y  if  \Jji{x,p)\^0, 

H 

0  otherwise.  (1) 


where  U  is  the  number  of  non-zero  entries  in  row  i  (i.e.  the 
number  of  components  for  which  component  i  influences  the 
change  of  state).  Note  that  an  arbitrary  state-dependent  matrix 
A(x)  could  be  transformed  to  a  stochastic  matrix  this  way. 
From  the  matrix  M  we  can  decompose  the  system  to  its 
horizontal-vertical  decomposition  (HVD)  [9].  The  decompo¬ 
sition  is  obtained  by  spectral  analysis  of  M.  This  is  the  same 
as  the  one  used  in  [3]  except  that  the  decomposition  itself  is 
related  to  Markov  chain  theory  and  different  vertical  levels 
are  identified  as  left  eigenvectors  of  M.  For  completeness 
and  to  set  the  decomposition  in  reader’s  mind  we  show  here 
a  figure  from  [9]  in  which  the  decomposition  is  graphically 
depicted.  The  graph  theoretic  decomposition  -  equivalent  to 
the  one  described  by  [3]  is  provided  in  the  next  section. 
These  decompositions  split  the  system  into  vertical  levels,  with 
components  on  each  level  affecting  change  only  in  the  ones 
above  them.  This  is  very  significant  for  uncertainty  propaga¬ 
tion.  Assume  there  is  parametric  uncertainty  in  parameter  p 
that  only  affects  a  state  at  the  lowest  level.  Since  dynamics  at 
the  lowest  level  is  not  affected  by  the  dynamics  at  the  levels 
above,  the  asymptotic  probability  density  at  the  lowest  level 
can  be  used  as  the  input  for  the  states  above.  In  this  way, 
the  computational  effort  for  probability  density  propagation  is 
substantially  reduced. 


A.  Decomposition  based  on  graph  theory 

We  now  present  an  alternative  method  to  obtain  the  HVD  of 
a  dynamical  system  using  graph  theory.  The  key  observation 
is  that  the  recurrent  states  obtained  recursively  using  the 
Markov  chain  theory  for  the  HVD  correspond  to  the  strongly 
connected  components  (SCC)  of  the  equation  graph.  Note  that 


Fig.  1.  Decomposition  described  in  the  theorem 


any  digraph  can  be  uniquely  decomposed  into  its  SCC  [5], 

[10]. 

Each  SCC  of  the  equation  graph  can  be  identified  as  a 
subsystem  (a  supemode)  and  the  condensed  graph  represent¬ 
ing  the  interactions  between  the  supemodes  is  an  acyclic 
digraph  (ADG)  [5].  The  HVD  can  be  obtained  in  a  straight 
forward  manner  from  this  condensed  graph.  The  subsystems 
corresponding  to  vertices  with  out-degree  zero  belong  to  level 
1  of  the  HVD.  These  vertices  can  be  deleted  and  the  process 
repeated  to  successively  identify  the  subsystems  in  the  next 
level. 

The  SCC  of  a  digraph  can  be  identified  efficiently  using  the 
depth  first  search  (DFS)  algorithm  [5],  [10].  In  particular,  the 
following  graph  algorithm  generates  the  SCC:  i)  form  a  DFS 
forest  of  the  given  digraph  and  post-order  the  vertices,  ii)  form 
a  DFS  forest  of  the  reverse  graph  of  the  original  graph  {i.e., 
with  all  the  edges  reversed  in  direction)  using  the  ordering 
obtained  in  step  i).  The  SCC  are  given  by  the  DFS  forest  in 
step  ii).  The  DFS  of  a  graph  can  be  carried  out  using  efficient 
algorithms  and  the  computational  effort  is  of  the  order  n  +  e 
where  n  is  the  number  of  vertices  and  e,  the  number  of  edges 
[5]. 

HI.  Decomposition  based  on  the  Jacobian  value 

In  this  section  we  investigate  the  use  of  the  Jacobian  matrix 
to  decompose  a  system  of  equations  into  a  collection  of  weakly 
coupled  systems,  each  of  which  consists  of  strongly  coupled 
equations.  We  note  that  the  Jacobian  matrix  is  generally  state- 
dependent,  hence  varies  with  time  as  the  states  evolve.  This 
variation  may  reflect  the  fact  that  in  different  regions  of  state 
space  the  strength  of  the  couplings  between  variables  change. 

Our  method  is  along  the  lines  of  [4]  where  transitions 
between  regions  in  state  space  for  a  dynamical  system  are 
described  by  a  Markov  chain  whose  underlying  graph  was 
employed  to  identify  the  number  of  almost  invariant  sets.  A 
symmetrized  adjacency  matrix  is  formed  from  the  stochastic 


matrix  of  the  Markov  chain  and  from  this  the  graph  Laplacian 
is  constructed.  The  eigenvalues  of  the  graph  Laplacian  are  used 
to  determine  the  number  and  location  of  the  almost  invariant 
sets  in  state  space. 

Here  we  propose  the  use  of  a  normalized  numerical  Jacobian 
matrix  to  replace  the  Markov  chain  probability  transition 
matrix.  Then  we  proceed  analogously  to  build  an  equation 
graph.  Eigenvalues  of  the  associated  graph  Laplacian  are  used 
to  identify  the  number  of  weakly  coupled  subsystems. 

Given  an  n  dimensional  system, x  =  f(x),  we  form  the 
matrix  W  (transpose  of  the  Jacobian)  with  entries  evaluated 
at  a  given  point  in  state  space 


dfj  (x) 

dxi 


The  entries  of  the  zth  row  are  normalized  by  their  row  sums 


n 


■&  =  £ 

3  =  1 


dfj  (g) 

dxi 


Normalizing  the  entries  by  their  row  sums,  we  get 


n 


3- 1 


dfj  0) 

dxi 


Let  the  normalized  matrix  be  P.  Clearly,  P  is  a  stochastic 
matrix,  and  as  such  can  be  thought  of  as  the  probability 
transition  matrix  for  a  Markov  chain.  In  our  method,  the  xi 
act  as  labels  for  the  graph  vertices.  The  P  matrix  measures 
the  ’effect’  that  Xi  has  on  each  Xj.  Next  we  form  the  diagonal 
matrix  M  =  diag{pi},  where  pi  is  the  zth  component  of 
the  stationary  probability  vector  for  P.  Finally,  we  form  the 
symmetric  matrix 


A  =  i  (MP  +  ( MP)T ) 

and  use  it  as  a  weighted  adjacency  matrix  for  our  equation 
graph,  with  vertices  labeled  by  the  states  Xi.  This  sym- 
metrization  corresponds  to  converting  a  directed  graph  into  an 
undirected  graph,  or  equivalently,  to  reversibilizing  a  Markov 
chain. 

In  order  to  determine  the  number  of  subsystems  the  system 
should  be  decomposed  into,  we  form  a  generalized  Laplacian 
matrix  (see  [6])  associated  with  the  graph  determined  by 
the  weighted  adjacency  matrix,  A.  To  form  the  generalized 
Laplacian,  we  first  compute  the  degree  of  each  vertex  in  the 
graph  determined  by  A  as 


deg(vj)  =  ^  Aij 
j= i 

where  each  vertex  Vi  corresponds  to  one  of  the  state 
variables  Xi.  Finally  we  form  the  generalized  Laplacian  matrix 

L  =  D~ ?  (D-A)D~^. 
where  D  =  diag  {deg  (u*)}  is  the  degree  matrix. 


Air  from 


Fig.  2.  Schematic  of  the  fuel  processing  system  of  a  fuel  cell  power  plant 
showing  various  reaction  stages 


Note  that  L  is  symmetric,  positive  semidefinite,  and  so  has 
no  negative  eigenvalues.  The  multiplicity  of  the  zero  eigen¬ 
value  gives  (see  [6])  the  number  of  connected  components 
of  the  graph  determined  from  A.  Numerically,  while  we  will 
always  find  one  zero  eigenvalue,  there  may  be  others  that  are 
close  to  zero.  In  [4],  the  multiplicity  of  the  zero  (and  near-zero) 
eigenvalues  is  used  as  an  indicator  of  the  number  of  almost 
invariant  sets  in  the  phase  space.  In  our  case,  the  connected 
components  of  the  graph  map  to  subsystems  of  equations. 
Consequently  we  use  the  multiplicity  of  the  zero  (and  near¬ 
zero)  eigenvalues  to  determine  the  number  of  weakly  coupled 
subsystems  into  which  our  original  system  decomposes.  An 
example  is  shown  in  the  next  section. 

As  mentioned  previously,  with  the  evolution  of  the  dy¬ 
namical  system,  the  state  variables  change  and  so  does  the 
evaluated  Jacobian  matrix.  This  means  that  the  decomposition 
may  be  different  in  different  regions  of  state  space.  The 
current  approach  is  to  mimic  the  modified  Newton  approach  to 
nonlinear  equation  solvers,  where  a  new  Jacobian  is  evaluated 
when  certain  convergence  criteria  are  not  met. 

IV.  Application  to  a  fuel  cell  system  model 

In  this  section,  we  demonstrate  the  two  decomposition 
algorithms  discussed  in  Sections  II  and  III  using  a  fuel  cell 
system  model.  The  fuel  cell  system  consists  of  a  fuel  processor 
that  combines  natural  gas  fuel  with  air  and  reforms  it  into 
hydrogen  rich  gas  using  a  series  of  reactors.  A  schematic  of 
the  fuel  processor  (without  the  fuel  cell  itself)  is  shown  in 
Figure  2. 

The  hydrogen  rich  gas  feeds  into  the  anode  channel  of  a 
PEM  (polymer  electrolyte  membrane)  fuel  cell  which  com¬ 
bines  hydrogen  from  the  anode  with  air  from  the  cathode  to 
produce  electricity,  heat  and  water.  The  model  we  consider 
has  highly  simplified  physics  and  was  developed  for  control 
studies  [11],  [12].  The  model  has  10  states  and  two  control 
inputs.  Figure  3  shows  the  lumped  volumes  in  the  model  and 
the  respective  state  variables  in  the  volumes. 

The  pressures  and  partial  pressures  in  the  volumes  are 
denoted  by  p.  The  superscripts  hex,hds,wrox,an  refer  to  the 
volume  (component)  name  and  the  subscripts  H2 ,  CH4  and 
air  refer  to  the  species.  The  temperature  of  the  CPOX 
reactor  is  denoted  by  Tcpox  and  the  speed  of  the  blower 
by  The  states  of  the  system  are  ordered  as  x  = 


Fig.  4.  Equation  graph  of  the  fuel  cell  system  model 


f rp  „ an  ^an  nhex  ,  ^hds  „ mix  mix  wrox  nwroxV 

[J-epox)  Ph?>  P  i^bloiP  i  Pch4'>  ^air  >Ph2  j  “  J  ' 

The  equation  graph  of  the  system  showing  the  who-affects- 
who  relation  is  shown  in  Figure  4.  The  Jacobian  of  the  system 
at  an  operating  point  of  interest  is  given  in  Eq  2 

Applying  the  Markov  chain  decomposition  method  of  Sec¬ 
tion  ??,  we  obtain  the  HVD  as  shown  in  Figure  5.  At  levels 
1,  2  and  4,  we  have  subsystems  with  a  single  state  whereas  at 
level  3,  we  have  a  subsystem  with  7  states. 

Using  the  graph  theoretic  algorithm  described  in  Section  II- 
A,  a  DFS  tree  for  the  equation  graph  can  be  constructed  with 
post-ordering  of  the  vertices  as  shown  in  Figure  6(a).  The 
DFS  forest  for  the  reverse  graph  using  the  ordering  from 
the  previous  step  is  shown  in  Figure  6(b).  This  yields  the 


Fig.  6.  (a)  DFS  tree  of  the  equation  graph  for  the  fuel  cell  system  with  post- 
ordering  of  the  vertices,  (b)  The  DFS  forest  of  the  reversed  graph  showing 
the  strongly  connected  components  of  the  original  graph 

strongly  connected  components  of  the  original  equation  graph. 
Notice  that  the  strongly  connected  components  are  precisely 
the  same  as  the  subsystems  at  various  levels  in  Figure  5  from 
the  Markov  chain  method. 

The  subsystem  at  level  3  in  Figure  5  contains  7  states 
and  is  a  strongly  connected  component.  The  Jacobian  value 
method  described  in  Section  III  can  be  applied  to  further 
decompose  this  subsystem  into  weakly  connected  subsystems. 
We  construct  the  stochastic  matrix  and  the  Laplacian  of  the 
associated  graph.  The  eigenvalues  of  the  Laplacian  are 

[  0  0.0296  0.1958  0.6518  0.8769  0.9912  1.1089  ] 


Level  1 


Level  2 


Level  3 


Level  <4 


Fig.  5.  Horizontal-vertical  decomposition  of  the  fuel  cell  system  model 


and  shown  in  Figure  7.  It  can  be  observed  that  the  second 
eigenvalue  is  relatively  closer  to  zero  than  the  rest  of  the 
eigenvalues.  This  indicates  that  the  7  state  system  can  be 
decomposed  into  two  weakly  coupled  subsystems.  The  two 
subsystems  can  be  identified  as  {TcpoX)pan}pWROX}  and 
{phex  7  phds  ?  pgj* ,  p7^  }  based  on  the  sign  of  the  components 
of  the  second  eigenvector.  Physically,  this  decomposition  re¬ 
lates  to  the  downstream  and  upstream  dynamics  of  the  CPOX 
reactor. 

V.  Conclusion/Summary 

Promise  of  graph  decomposition  methods  for  enabling 
efficient  uncertainty  propagation  in  complex,  networks  of 
dynamical  systems. 
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Spectral  balance:  a  frequency  domain  framework  for  analysis  of  nonlinear  dynamical  systems 


Andrzej  Banaszuk  Prashant  G.  Mehta 


Abstract — A  frequency-domain  framework  for  analysis, 
computations,  and  uncertainty  propagation  in  nonlinear  sys¬ 
tems  driven  by  broad-band  disturbances  is  introduced  and 
illustrated  in  a  simple  example  of  a  nonlinear  system  that 
exhibits  noise-induced  transitions  between  two  stable  equilib¬ 
ria.  The  spectral  balance  framework  generalizes  the  standard 
harmonic  and  Gaussian  signal  balance  in  feedback  systems. 
The  application  example  presented  is  a  scalar  model  with 
cubic  nonlinearity  after  pitchfork  bifurcation  driven  by  a 
broad-band  disturbance.  An  approximate  and  iterative  spec¬ 
tral  balance  (including  determination  of  equilibria)  is  solved. 
The  solution  of  the  approximate  spectral  balance  is  used  to 
reformulate  the  original  model  using  a  loop  transformation 
so  that  an  iterative  procedure  for  finding  the  spectrum  of  the 
output  converges  to  the  true  spectrum  of  the  solution. 

I.  Introduction 

Many  industrial  flows  involve  complex  interactions  of 
acoustic  waves,  vorticity,  fuel  transport,  and  chemical  re¬ 
actions.  The  control  objective  often  is  to  create  beneficial 
non-equilibrium  dynamics  with  control.  Examples  include 
control  of  flow  separation  and  mixing  enhancement.  In 
this  paper  we  introduce  a  frequency  domain  framework 
for  analysis  and  non-equilibrium  control  design  for  a  large 
class  of  models  of  physical  phenomena  involving  mul¬ 
tiple  oscillatory  modes  coupled  through  nonlinear  terms, 
transport  delay,  and  driven  by  broad-band  disturbances. 
While  motivated  by  specific  problems  arising  in  military 
aeroengines,  the  methods  will  be  applicable  to  large  class  of 
distributed  dynamical  systems  involving  oscillatory  dynam¬ 
ics  with  nonlinear  cross-coupling,  saturated  nonlinearities, 
transport  delay,  and  broad-band  disturbances. 

The  spectral  balance  framework  that  we  propose  gener¬ 
alizes  the  standard  harmonic  balance  and  Gaussian  signal 
balance  in  feedback  systems  [4],  [2].  The  framework  is 
introduced  and  illustrated  in  an  example  of  a  nonlinear 
system  that  exhibits  noise-induced  transitions  between  two 
stable  equilibria.  The  example  presented  is  a  scalar  model 
with  cubic  nonlinearity  after  pitchfork  bifurcation  driven 
by  a  broad-band  disturbance.  An  approximate  and  iterative 
spectral  balance  of  the  constant  and  broad-band  signals  (in¬ 
cluding  determination  of  equilibria)  is  solved.  The  solution 
of  this  approximate  spectral  balance  is  used  to  reformulate 
the  original  model  using  a  loop  transformation  so  that  an 
iterative  procedure  for  finding  the  spectrum  of  the  output 
converges  to  the  true  spectrum. 
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II.  Spectral  Balance 

Consider  a  model  of  a  lightly  damped  stable  linear  system 
with  transfer  function  Go(joj),  in  a  feedback  loop  with  a 
static  nonlinearity  /(•),  subject  to  a  driving  disturbance  n(t) 
with  the  Fourier  transform  N(ju).  An  uncertainty  in  the 
model  is  represented  by  an  (in  general  nonlinear)  operator 
A(-)  in  a  feedback  loop  around  the  nominal  model  The 


Nominal  Model 


Fig.  1.  The  model  structure 
model  equations  are 

X(juj)  =  G0(ju>)(N(jw)  -  A(juj,X(juj))  -  Y(ju>))  (1) 

y(t)  =  f(x{t))  (2) 

where,  X(')  =  !Fx{-),  Y(-)  =  Xy{-),  and  N(-)  =  JYi(-), 
are  the  Fourier  transforms  of  the  corresponding  temporal 
signals.  We  assume  that  nonlinear  mapping  /(•)  is  Lipshitz 
on  each  bounded  set.  The  equation  (2)  can  be  represented 
in  the  frequency  domain  as 

Y(ju)  =  f(X(ju>))  :=  Tf{F~lX{ju)).  (3) 

Now,  the  feedback  system  (l)-(2)  can  be  represented  as 

X(juj)  =  GoUu){N(ju)  -  f(X(ju)  -  A(ju)X(ju)). 

(4) 

Note  that  for  the  linear  case  f(x)  =  0,  A (ju,X(juj))  = 
A (jLo)X(joj)  the  mapping  of  the  uncertainty  A^u;)  to  the 
output  of  the  system  is  given  by  the  formula 

X(joj)  =  (I  +  GoU^Aiju^GoU^WM.  (5) 

involving  the  sensitivity  function  (I  +  Gotjcj) AO'u;))-1. 
Note  that  the  frequency  domain  representation  greatly 


simplifies  the  uncertainty  propagation  analysis. 

1.  Uncertainty  propagation.  The  sensitivity  function 
(I  +  G0(j^)A(jo;))“1  allows  to  explicitly  map  the 
probability  distribution  of  the  uncertain  parameters 
contributing  to  A (ju>)  to  the  probability  distribution  of  the 
output  Y(ju). 

2.  The  frequency  domain  representation  greatly  accelerates 
computation  of  this  mapping.  Note  that  only  the  algebraic 
calculations  need  to  be  performed  in  evaluating  the  formula 
(5).  In  contrast,  a  time  domain  counterpart  of  the  (5)  would 
require  evaluation  of  the  convolution  integrals  over  long 
period  of  time.  Moreover,  in  the  time  domain  formulation 
one  needs  to  wait  for  transients  to  subside,  which  is 
the  issue  when  dealing  with  lightly  damped  dynamics. 
There  are  additional  benefits  of  the  frequency  domain 
representation  in  the  case  when  Go(ju)  contains  time 
delays. 

3.  Tools  from  the  robust  linear  control  theory  allow 
to  handle  dynamic  uncertainty  in  case  when  only  the 
bounds  on  the  uncertain  operator  are  known  [1]. 

Appart  from  application  to  the  uncertainty  propagation, 
the  frequency  domain  formulation  allows  to  study  funda¬ 
mental  limitations  of  achievable  control  performance  using 
methods  of  the  complex  analysis. 

The  spectral  balance  approach  is  a  frequency  domain 
framework  for  the  uncertainty  analysis  of  nonlinear  systems 
that  includes  the  case  of  non-equlibrium  bounded  dynamics 
that  retains  the  advantages  of  the  linear  sensitivity  function 
framework:  explicit  formulas  mapping  uncertainty  to  the 
output  and  the  speed  of  computation.  We  will  begin  with 
the  particular  case  of  system  in  Figure  1  using  the  fixed 
point  formulation  (4).  To  introduce  the  spectral  balance 
framework  we  will  consider  the  case  without  uncertainty 
shown  in  Figure  2  with  the  corresponding  fixed  point 


Fig.  2.  The  model  structure 


formulation  of  the  spectral  balance  equation  given  by  (6) 

X(joj)  =  G0(joj)(N(juj)  -  f(X(ju)).  (6) 

Note  that  the  spectral  balance  framework  generalizes  the 
standard  harmonic  balance  (where  the  input  signal  n(-) 


is  periodic,  or  when  the  dynamics  has  limit  cycles)  and 
Gaussian  signal  balance  (where  the  input  signal  n(-)  is  a 
Gaussian  broad  band  signal)  in  feedback  systems  [4],  [2], 

We  assume  that  the  dynamics  of  (2)  is  globally  bounded 
and  that  there  is  an  attractor.  Eventually  we  intend  to 
introduce  a  spectral  balance  framework  for  the  class  of 
bounded  power  signals  on  an  infinite  time  interval.  In  this 
paper  we  restrict  the  attention  to  the  space  of  L2  signals  on 
the  interval  [0,  T],  where  T  is  large  relative  to  the  slowest 
time  scale  in  the  system.  The  induced  operator  norms  are 
the  Hr^  norms. 

A  sufficient  condition  for  existence  of  a  unique  solution 
of  the  spectral  balance  equation  (6)  is 

||G0(^)(/(X2(ja;))-/(X1(^)))||  <  WX^-X^W- 

(7) 

For  all  Xi(juj)  in  I/2[0,  T}.  Note  that  in  this  case  a  unique 
solution  to  (6)  exists  (by  applying  the  Banach  Contraction 
Mapping  Theorem  [3]).  Moreover,  the  approximate  solu¬ 
tion  of  the  spectral  balance  equation  can  be  obtained  by 
successive  approximations  using  the  formula 

Xi+i(ju)  =  G0(juj)(N(ju>)  -  f(Xi(ju)).  (8) 
with  an  arbitrary  initial  condition. 

Note  that  if  the  condition  (7)  is  satisfied  for  all  Xi(ju)  in 
a  closed  set  B  in  L2[0,  T\  that  is  invariant  for  the  mapping 
Go(ju)(N (ju)  -/(*)),  one  can  approach  a  solution  of  the 
spectral  balance  equation  (6)  in  B  using  (8)  with  Xo(jcj)  £ 
B . 


III.  LOOP  TRANSFORMATION 

A  sufficient  condition  for  (7)  is  the  small  gain  condition 
for  the  feedback  loop  in  (2).  However,  even  if  the  condition 
(7)  is  not  satisfied,  which  would  be  the  case  if  the  loop  gain 
is  large,  one  can  attempt  to  enforce  the  condition  (7)  for  an 
equivalent  feedback  system  to  (2)  by  a  loop  transformation. 

An  example  of  a  linear  loop  transformation  is  shown  in 
Figure  3.  Here  H(ju)  is  an  arbitrary  stable  linear  operator, 

G\ (ju)  :=  (I  +  H{ju)GQ{ju))-'G0(joj).  (9) 

and 

h(X(jw))  :=  /(*(,«))  -  H(ju)X(juj).  (10) 

The  spectral  balance  condition  for  the  system  in  Figure  3 
is 

X(ju)  =  Gi(ju){N{ju)  -  fi(X(M).  (11) 

Note  that  a  sufficient  condition  for  the  contraction  condition 
for  the  transformed  spectral  balance  (11)  is 

I  l^i  0w)  (A  (-^2  (Ju))~  h{Xi  (ju) ) )  1 1  <  wxzijuyx^Mi 

(12) 

If  the  nonlinear  part  of  the  loop  in  Figure  2  has  a  stabilizing 
effect,  the  role  of  the  operator  H(joj)  is  to  reduce  the 
gain  of  the  nonlinear  part  of  the  loop  and  increase  the 


Fig.  3.  The  loop  transformation  to  enforce  loop  contractivity 


contractivity  of  the  linear  part  of  the  loop.  More  precisely, 
the  gain  of  the  nonlinear  operator  f(X(ju)))  is  reduced  by 
subtraction  of  a  linear  approximation  of  and  the 

approximate  linear  operator  is  incorporated  in  the  linear  part 
of  the  modified  loop.  Thus,  a  good  choice  of  H(juj)  is  the 
one  that  minimizes  \\fi(X(juj))\\  for  X(ju)  representing 
the  solution  of  the  spectral  balance  equation  (11). 

Of  course,  the  minimization  of  \\fi(X(juj))\\  requires 
knowledge  of  X{ju)  itself,  which  is  exactly  the  solution  of 
the  spectral  balance  equation  that  we  seek.  Note  that  we  are 
interested  in  the  case  of  system  (2)  having  non-equilibrium 
attractors  or  subject  to  a  large  driving  disturbance,  so  that 
an  approximation  of  nonlinear  operator  f(X(ju))  by  its 
linearization  at  ( X(ju )  —  0  is  not  suitable. 

The  key  idea  introduced  in  this  paper  is  to  proceed  in 
the  following  three  steps: 

1.  Find  an  approximate  solution  Xappr(ju)  close  to 
X{ju).  For  this,  the  describing  function  techniques,  both 
for  harmonic  and  random  Gaussian  signals,  can  be  utilized. 
In  fact,  as  we  will  show  in  the  next  section,  it  may  be  not 
necessary  to  find  an  approximate  solution  Xappr{juj ),  but 
only  few  parameters  describing  such  a  solution,  like  its 
time  average  and  the  average  power. 

2.  Utilize  the  knowledge  of  Xappr{ju )  (or  the  parameters 
describing  it)  to  find  a  linear  transformation  H(jui)  that 
minimizes  (or  at  least  reduces)  \\fi(Xappr(juj))\\  = 
\\f(Xappr(joj))  -  H(joj)Xappr(ju) ||. 

3.  Use  H(juj)  to  define  the  loop  transformation  (9)- 

(10).  If  the  contraction  condition  (12)  is  satisfied  on  a 
closed  set  B  in  L2[0,T]  that  is  invariant  for  the  mapping 
Gi(ju)(N(ju)  -  one  can  approach  a  solution  of 

the  spectral  balance  equation  (6)  in  B  using  the  iterative 


process  defined  by 

*i+i (joj)  =  G ! -  h (Xi(ju))  (13) 

starting  with  an  arbitrary  Xq  (juS)  6  B. 

At  present  it  is  not  clear  under  what  general  conditions 
the  procedure  described  above  will  result  in  finding  solu¬ 
tions  to  the  spectral  balance  equations.  In  the  next  section 
we  will  show  one  example  of  a  system  with  nontrivial 
dynamics,  for  which  that  the  procedure  yields  the  desired 
result. 

IV.  Example 
Consider  the  equation 

x(t)  +  ax(t)  +  bx3(t)  =  n(t).  (14) 

Here  we  assume  that  a  <  0,  b  >  0,  and  the  input  signal 
n(-)  has  zero  mean  and  flat  spectrum  |jV(jiu;)|  =  cr*  for 
all  u).  In  the  sequel  we  will  refer  to  the  input  signal  n(-) 
as  noise ,  even  though  we  emphasize  that  in  this  paper  we 
only  consider  the  deterministic  case.  Note  that  for  a  =  0 
the  system  (14)  undergoes  a  pitchfork  bifurcation  and  the 
equilibrium  x  =  0  becomes  unstable  for  all  a  >  0.  Two 
locally  stable  equilibria  occur  at  x  =  and  x  = 

jsjote  that  for  a  small  value  of  of  the  solution 
x(t)  will  be  close  to  one  of  the  stable  equilibria.  For  some 
higher  value  of  of  the  solution  x(t)  will  be  transitioning 
from  a  neighborhood  of  the  one  of  the  stable  equilibria 
to  the  other,  as  shown  in  Figure  4.  Note  that  the  spectral 


Typical  solution 


Time  trace,  a=-1 ,  b=1 ,  noise  power=0.99472 


Time,  sec 


Fig.  4.  Typical  solution  of  system  (14):  noise-induced  transitions  between 
two  stable  equilibria 

balance  equation  for  (14)  is 

=  (N(j<s)-f(X(ju)),  (is) 

jio  +  a 


where  f(x)  :=  bx3.  Note  that,  since  a  <  0,  the  linear 
operator  is  unstable,  and  the  contraction  condition 

(7)  is  not  satisfied.  However,  the  nonlinear  operator  f{x) 
has  a  stabilizing  effect,  so  that  we  can  attempt  to  transform 
the  loop  to  an  equivalent  one,  for  which  the  contraction 
condition  is  satisfied,  as  described  in  Section  III. 

Let  x  :=  J  Jq  x(t)dt  denote  the  time  average  of  x(t) 
and  let  xf(t)  :=  x(t)  -  x  denote  the  deviation  of  x(t) 
from  its  average  value.  Let  <r2  :=  ^  JQT  x'(t)2dt  denote 
the  mean  power  of  x'(t).  In  what  follows  we  will  compute 
approximate  value  of  x  and  cr2  by  solving  approximate 
spectral  balance  equations.  By  taking  the  time  average  of 
(16)  and  using  the  fact  that  the  average  of  (x  +  x'(t))3  is 
x3  +  3  bal  we  obtain 


Now,  solving  (16)  and  (23)  we  obtain 


>  46 


x2  =  -f  -  3 a2 


a2  =  Hfi  +  ,/i  + 

°x  2b  v1  ^  V  ^  a2  h 

x  =  0. 


(24) 

(25) 

(26) 

(27) 

(28) 
(29) 


Figure  5  graphically  represents  the  solution  to  the  approx¬ 
imate  spectral  balance  equations  as  function  of  the  input 
power  <7;  for  a  =  —1,  b  =  1.  The  solutions  for  x  and 


x(a  +  3ba2  +  bx2)  —  n  =  0.  (16) 

Subtracting  (16)  from  (14)  and  re-arranging  terms  yields 
x'(t)  +  (a  +  h)x'(t)  +  fi(x'(t))  =  n(t),  (17) 

where 

h  :=  ba\  +  (18) 


/i(x'(t))  :=  b(x'2  -  a2x){x'  +  3x).  (19) 

To  find  approximate  values  of  x  and  a\  we  will  neglect  the 
term  h(x>  ax,x'(t))  in  (17)  and  solve  the  equation 

x'(t)  +  (a  +  bal  +  3  bx2)xf(t)  —  n(t).  (20) 


Note  that  for  fixed  values  of  x  and  a\  (20)  is  a  linear 
equation  that  can  be  solved  in  the  frequency  domain  as 


X'(jcj)  = 


N(ju) 

ju  +  a  -f  bal  + 


(21) 


For  a  moment  we  assume  that  the  values  of  x  and  a\  are 
such  that  a  +  bal  +  >  0,  so  that  the  transfer  function 

ju+a+ba2  +36P  us  sta^e-  This  assumption  will  be  verified 
after  x  and  a\.  are  calculated. 

Now  we  obtain  the  closure  equation  for  a\  by  integrating 
the  absolute  values  of  both  sides  of  (21)  over  all  frequencies 


*  2tt  J_c 


ju  +  a  +  bal  +  3bx2 


2  duo. 


(22) 


The  equations  (16)  and  (22)  form  an  approximate  spectral 
balance  for  the  system  (14.  The  integral  in  (22)  can  be 
analytically  evaluated  so  that  we  can  write  the  following 
equation 


2 (a  -b  bal  +  3frr2)  * 


(23) 


Approximate  spectral  balance 


Fig.  5.  Solutions  to  approximate  spectral  balance  equations  as  function 
of  the  input  power 

al  have  a  natural  intepretation.  For  there  are 

two  values  of  the  time  average  x  close  to  the  no-noise 
equilibria  that  can  be  attained.  The  value  of  ax  (that  could 
be  interpreted  as  standard  deviation  of  x(t))  is  small,  so 
that  the  solution  x(t)  stays  close  to  an  equilibrium  solution 
and  does  not  transition  to  the  neighborhood  of  the  other 
equilibrium.  Above  the  critical  value  of  the  input  power 
a{  =  the  solutions  x(t)  deviate  from  the  stable  equilibria 
far  enough  to  transition  between  the  neighborhoods  of  the 
both  equilibria.  Since  the  transitions  back  and  forth  can 
occur,  x  =  0  becomes  the  mean  and  the  standard  deviation 
ax  is  close  to  the  distance  from  the  new  mean  x  —  0  to 
the  value  where  the  solution  x(t)  spends  most  of  the  time: 
close  to  the  no-noise  equilibria  of  (14). 

We  will  now  use  the  values  of  x  and  a2x  given  by  (5)  to 
perform  a  loop  transformation  as  described  in  Section  III. 


More  precisely,  we  will  solve  the  perturbation  equation  (17) 
in  the  frequency  domain  using  the  fixed  point  formulation 

x'U»)  =  — -  A (30) 

ju)  +  a  H-  fi 

It  can  be  easily  verified  that  a-f/i  >  0  for  x  and  given  by 
(29).  Analytic  verification  of  the  contraction  condition  for 
the  operator  j~^a+h  -  /i(-)  is  difficult.  Therefore 

we  will  assume  that  the  contraction  condition  is  satisfied 
and  proceed  with  an  iterative  solution  to  (30)  using 

=  4l,.\.h(NV»)  -  fi(XW)  (31) 

juj  -H  a  +  h 

with  X'o(juj)  —  0  and  verify  the  contraction  condition 
numerically.  To  illustrate  and  verify  this  procedure,  a  nu¬ 
merical  solution  of  (14)  for  a  —  -1,  b  ~  1,  and  cr;  >  — 
was  obtained.  The  spectrum  N(ju)  of  the  noise  from 
the  time  domain  simulation  was  saved  and  used  in  the 
formula  (31).  Figure  6  shows  an  excellent  agreement  of 
the  spectrum  X'(jui)  from  the  time  domain  simulation 
and  the  spectrum  X'iq  (jui)  from  the  iterative  procedure 
(31)  after  10  iterations.  Figure  7  shows  comparison  of  the 


Simulation  and  approximation 
after  10  iterations:  time  traces 


Time  traces,  a=-1,  b=1,  noise  power=3.9789 


Fig.  7.  Solution  to  iterative  spectral  balance  equations:  time  traces 


Simulation  and  approximation 
after  10  iterations:  spectra  (fft) 


Spectra,  a=-1,  b=1,  noise  power=3.9789 
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Fig.  6.  Solution  to  iterative  spectral  balance  equations:  spectra 

time  traces  of  the  solutions  of  (14)  obtained  by  the  time 
domain  simulation  and  by  the  iterative  spectral  balance  and 
the  inverse  Fourier  transform.  Figure  8  shows  decay  of 
the  power  of  the  approximation  error  X'(jui)  -  X'i(juj) 
normalized  by  the  power  of  X'(jcj)  as  a  function  of 
,  Figure  9  shows  the  contraction  rate 
as  a  function  of  iteration  step  i.  This 
at  the  rate  of  about  0.8  was  indeed 
achieved  by  the  loop  transofrmation  involving  solution  of 
the  approximate  spectral  balance. 


iteration  step  i .  Finally. 

\\X'i+1{j»)-X'i{ju)}\ 
verifies  the  contraction 


210  FFT  points 


Relative  error  of  approximation,  a=-1 .  b=1 ,  noise  power=3.9789 


Number  of  iterations 


Fig.  8.  Solution  to  iterative  spectral  balance  equations:  relative  approxi¬ 
mation  error  for  210  FFT  points 

In  the  case  that  were  examined  obtaining  an  approximate 
solution  of  (14)  using  the  formula  (31)  was  orders  of 
magnitude  faster  that  the  time  domain  simulations  using 
Simulink. 

V.  Conclusion 

A  frequency-domain  framework  for  analysis,  computa¬ 
tions,  and  uncertainty  propagation  in  nonlinear  systems 
driven  by  broad-band  disturbances  was  introduced  and 
illustrated  in  a  simple  example  of  a  nonlinear  system 
that  exhibits  noise-induced  transitions  between  two  stable 
equilibria.  The  spectral  balance  framework  generalizes  the 
standard  harmonic  and  Gaussian  signal  balance  in  feedback 


Contraction  rate,  a=-1 ,  b=1 ,  noise  power=3.9789 
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Fig.  9.  Solution  to  iterative  spectral  balance  equations:  contraction  rate 


systems.  The  application  example  presented  is  a  scalar 
model  with  cubic  nonlinearity  after  pitchfork  bifurcation 
driven  by  a  broad-band  disturbance.  An  approximate  and 
iterative  spectral  balance  (including  determination  of  equi¬ 
libria)  is  solved.  The  solution  of  the  approximate  spectral 
balance  is  used  to  reformulate  the  original  model  using 
a  loop  transformation  so  that  an  iterative  procedure  for 
finding  the  spectrum  of  the  output  converges  to  the  true 
spectrum  of  the  solution.  The  future  work  will  involve  more 
carefull  study  of  the  function  spaces  suitable  for  the  spectral 
balance  formulation  and  obtaining  some  analytic  sufficient 
conditions  for  the  contraction. 
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Abstract — Efficient  propagation  of  probability  distributions  in 
feedback  systems  is  central  to  a  decomposition  based  approach 
for  uncertainty  propagation  in  complex,  interconnected  systems. 
In  this  note,  we  propose  an  iterative  method  for  static  feedback 
systems  to  obtain  the  probability  density  of  the  output  from  that 
of  the  input.  We  prove  the  convergence  of  the  proposed  method 
under  the  assumption  that  the  loop  operator  is  contractive.  The 
method  is  illustrated  with  an  example.  It  is  shown,  based  on 
the  results  from  the  theory  iterated  random  functions,  that 
the  method  extends  to  the  case  when  additional  parametric 
uncertainty  is  present  within  the  loop. 


I.  Introduction 


Robust  design  of  complex,  interconnected  systems  requires 
efficient  computational  methods  for  uncertainty  propagation. 
Knowledge  of  the  probability  distribution  function  (pdf)  of 
key  output  variables  derived  from  the  knowledge  of  input  vari¬ 
ables  enables  better  decision  making  during  design.  Existing 
methods  for  propagation  of  uncertainty  such  as  Monte-Carlo, 
polynomial  chaos  and  stochastic  surface  response  methods 
have  scale  poorly  with  problem  size  and  are  computationally 
very  demanding  for  complex  systems  [1]. 

In  order  to  overcome  the  barrier  of  computational  effort, 
a  new  approach  based  on  decomposition  of  complex  systems 
has  been  emerging  as  a  promising  direction.  In  this  approach, 
a  complex  system  is  first  decomposed  using  graph  theoretic 
methods  into  subsystems  connected  in  series,  parallel  or  feed¬ 
back  [2].  Then,  a  block-by-block  propagation  is  performed 
accounting  for  the  dependencies  among  variables  [3].  Such 
a  method  exploits  the  underlying  hierarchical  structure  of  a 
complex  system  and  provides  the  flexibility  to  use  different 
methods  for  different  subsystems  of  the  original  system. 

Feedback  loops  at  the  system  level,  or  encompassing  a 
large  number  of  subsystems  pose  a  significant  risk  to  the 
block-by-block  framework  since  the  structural  decomposition 
method  described  in  [2]  identify  the  feedback  loops  as  a  single 
subsystem.  The  decomposition  method  based  on  the  Jacobian 
value  and  the  associated  graph  Laplacian,  also  described  in 
[2],  can  identify  weak  feedback  connections  only.  However, 
moderate  to  strong  feedback  ( e.g due  to  a  controller)  is 
common  in  many  engineering  systems.  The  ability  to  prop¬ 
agate  uncertainty  through  such  feedback  loops  block-by-block 
without  having  to  solve  the  closed  loop  system  in  entirety  is 
very  attractive.  Thus,  there  is  a  need  for  alternative  approaches 
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to  speeding  up  the  computation  for  uncertainty  propagation  in 
feedback  systems. 

In  this  note,  we  seek  an  iterative  method  to  compute  the 
pdf  of  the  output  of  a  feedback  loop  from  the  pdf  of  the 
input.  The  feedback  system  under  consideration  is  shown  in 
Figure  1.  The  loop  has  input  u,  the  internal  signal  x  and  output 
y  —  Gx  where  G  is  the  loop  operator  (possibly  nonlinear). 
In  Section  II,  we  present  the  problem  setup  and  formulate 
the  iteration  equations  by  abstracting  a  computational  scheme 
conceived  and  implemented  by  Kalmar-Nagy  and  Huzmezan 
[3].  We  prove  the  point- wise  convergence  of  the  iteration  to 
the  true  closed  loop  pdf  under  the  assumption  that  the  loop 
operator  is  contractive. 

In  Section  V,  we  consider  the  case  when  there  is  additional 
parametric  uncertainty  in  the  loop  operator.  We  show  that  the 
problem  can  be  cast  as  a  system  of  iterated  random  functions. 
Based  on  the  results  in  [4],  [5],  we  claim  that  the  iteration 
converges  to  the  solution  of  the  closed  loop  pdf  under  the 
assumption  that  the  loop  operator  is  contractive  on  an  average. 
The  proof  of  this  claim  and  extension  to  dynamic  systems  will 
be  considered  in  future  work. 

II.  Problem  formulation  and  iteration  scheme 

Consider  the  feedback  system  shown  in  Figure  1.  The  input 
u  belongs  to  a  Banach  space  X  (e.g.,  Rn)  and  is  randomly 
distributed  with  probability  measure  (iu.  The  assumption  of 
a  normed  linear  space  is  made  to  simplify  the  notions  of 
derivative  and  contraction  of  the  loop  operator,  but  it  may 
be  possible  to  relax  it  to  a  metric  space. 

Assuming  jiu  is  absolutely  continuous  with  respect  to 
the  Lebesgue  measure,  a  probability  density  function  (pdf), 
pu(u)  can  be  defined  for  u.  The  loop  operator  G  is  possibly 
nonlinear  and  maps  X  to  itself.  Furthermore,  we  assume  that 
G  is  contractive  in  the  sense  that  G  shrinks  distances  in  the 


associated  norm.  We  also  assume  that  G  is  one-to-one  and 
the  derivative  of  G ,  dG  exists  and  satisfies  ||dG||  <  1  which 
implies  that  G  is  uniformly  Lipschitz  with  Lipschitz  constant 
L  <  1.  Notice  that  the  assumption  of  a  single  loop  operator  G 
is  not  restrictive  since  G  can  be  thought  of  as  a  composition 
of  several  operators. 

The  problem  is  to  compute  the  pdf,  py{y)  of  the  output  y 
in  an  iterative  fashion,  without  directly  employing  the  closed 
loop  solution  y  =  G(J  -  G)~lu.  For  practical  reasons,  we 
restrict  our  treatment  to  finite  approximations  of  the  pdf’s  (i.e., 
discretized  versions).  This  allows  a  gridding  (box  covering) 
of  the  space  X  and  a  finite  matrix  representation  for  dG. 
Recently,  efficient  computational  methods  using  such  box  cov¬ 
erings,  known  as  set-oriented  numerical  methods,  have  been 
developed  and  extensively  applied  for  exploring  the  almost- 
invariant  sets  and  invariant  measures  of  dynamical  systems 
[6],  [7]. 

We  formulate  the  iterative  scheme  for  obtaining  px(x). 
Since  y  =  G(x),  py(y)  can  be  readily  obtained  from  the  pdf  of 
x ,  px(x)  using  the  Perron-Frobenius  theorem  [8,  chapter  17]: 


Py(v)  = 


Px{G~1y) 

\dG\x—G~ly 


Here  the  subscript  on  dG  indicates  where  the  derivative  is 
to  be  evaluated  and  |  •  |  denotes  the  absolute  value  of  the 
determinant.  Note  that  due  to  the  one-to-one  assumption  on 
G,  there  is  only  one  pre-image  for  y.  Otherwise,  the  right  hand 
side  would  involve  the  sum  over  all  pre-images  [8]. 

Using  the  Perron-Frobenius  theorem  on  the  closed  loop 
relation  x  =  (/  -  G)"1^,  the  pdf  of  x  is  given  by 


ptl(x)=Pu((I-G)x)\I-dG\  (1) 


We  define  below  a  sequence  iteration  for  px(x )  that  converges 
to  the  closed  loop  expression  in  Eq  1 .  The  proof  of  conver¬ 
gence  is  given  in  Section  III. 

We  first  describe  the  iterative  scheme  proposed  by  Kalmar- 
Nagy  and  Huzmezan  for  a  simple  scenario  [3]  .  Assume  that 
u,x,y  are  scalars  and  G  is  a  static,  one-to-one  map.  The 
proposed  iteration  scheme  grids  u  into  intervals,  say  of  length 
5u,  which  leads  to  a  discretization  of  the  pdf  as  a  histogram. 
The  first  iterate,  px{x)  is  taken  as  pu{x),  with  5xQ  =  6u.  Then, 
py(y)  is  given  by  p* y *  with  Sy°  =  \dG\Sx°.  The  second 
iterate,  plx{x)  is  obtained  by  making  use  of  the  fact  that  the 
area  under  the  curve  is  conserved,  ie,  pu(u)5u  —  p1x(x)5xl. 
The  interval  dx1  is  taken  as  the  sum  Su  +  6y°.  The  process 
is  repeated  until  px(x)  converges.  The  iterations  are  described 
by 


xk+1  =  u  +  Gxk 

(2) 

(5a;fc+1  =  6u  +  Syk 

(3) 

5yk  =  \dG\Sxk 

(4) 

,k+ iM  _  Pu{u)6u 
x  VX)  6xk+ 1  • 

(5) 

When  G  is  not  one-to-one,  the  ’’overlaps  of  measure”  in  the 
x  space  should  be  taken  care  by  adding  them  properly. 


We  now  present  a  generalization  of  the  above  iteration 
scheme  and  prove  its  convergence  under  the  assumption  that 
G  is  contractive  and  one-to-one.  Let  dG  denote  the  derivative 
of  G,  M  —  dxjdu  and  |  •  I  denote  the  absolute  value  of  the 
determinant. 

The  main  idea  of  the  iterative  method  is  to  define  a  discrete 
dynamical  system  from  the  given  static  feedback  system 
such  that  the  stationary  distribution  of  the  Perron-Frobenius 
equation  for  the  discrete  dynamical  system  is  the  same  as  the 
desired  closed  loop  distribution  for  the  static  feedback  system. 
The  key  steps  of  method  are  as  follows: 

1)  Partition  the  u  space  into  a  collection  of  boxes  Bi  with 
center  points  Ui 

2)  Assign  the  measure  pi  to  the  box  Bi  from  the  pdf  of  u , 

/.e,,  pi  -  pu{Bi) 

3)  Compute  the  image  of  each  Bi  iteratively  under  the 
assumption  that  the  image  set  can  be  characterized  by  x 
and  M  =  dxjdu .  This  holds  for  sufficiently  small  Bi. 

4)  Assign  the  measure  pi  to  the  (converged)  image  set  in 
x  space. 

The  underlying  assumption  is  that  if  the  size  of  Bi  is  suffi¬ 
ciently  small,  the  image  set  of  Bi  is  also  a  box  completely 
characterized  by  its  center  point,  x  and  the  linear  operator, 
M  =  dxjdu  that  captures  the  volume  mapping  from  the  u 
space  to  the  x  space.  The  probability  density  is  assumed  to  be 
uniform  with  in  each  box  for  all  variables. 

The  following  two  maps  are  employed  for  the  iteration: 

xk+l  =  u  +  Gxk  (6) 

Mk+l  =  I  +  (dG)Mk.  (7) 


The  initial  conditions  are  x°  =  u  and  M°  =  I.  The  iteration 
of  the  center  points  of  the  boxes  is  captured  by  Eq  6  and  the 
expansion  of  volume  elements  is  captured  by  Eq  7  which  is 
the  iteration  for  M  =  dxjdu. 

The  pdf  of  x  is  successively  approximated  as 


Px{xk)  = 


Puju) 
\Mk\  ' 


by  P(yk)  =  sSSr1- 


III.  Proof  of  Convergence 


Notice  that 

k 

Xk  =  (53  G>  =  {I  +  G  +  G2  +  ...  +  Gk)u. 

i=0 

If  G  is  linear,  i.e.9  if  dG  is  independent  of  x9  a  similar 
expression  can  be  obtained  for  Mk : 

k 

Mk  =  J2dGi  =  1  +  dG  +  d°2  +  ■  ■  ■  +  d°k ■ 

t=0 

However,  we  do  not  assume  linearity  of  G. 

Since  G  is  contractive,  by  the  contraction  mapping  theorem, 
the  map  in  Eq  6  converges  to  x*  =  (I  -  G)~lu.  Since  dG  is 


Fig.  2.  Convergence  of  the  iterates  of  pdf  of  x  for  the  feedback  loop  with 

y/x 

also  a  contraction,  the  map  in  Eq  7  converges  to  M*  =  (I  - 
dG)~l  with  dG  evaluated  at  x*.  Thus  the  pdf  of  x  converges 
to 

Px{x*)  =  =  Pu((I  -  G)x)\I  -  dG\ 

which  is  the  correct  closed  loop  solution  given  in  Eq  1 . 

Thus,  we  have  shown  point-wise  convergence  of  the  es¬ 
timated  pdf  of  x  to  the  closed  loop  pdf,  p^l{x)  at  selected 
grid  points.  A  well-known  result  due  to  Li  (see  [9]  and  [8, 
chapter  1 7])  guarantees  the  convergence  of  finite  approxima¬ 
tions  of  pdf’s  to  the  continuous  ones  as  the  grid  size  becomes 
sufficiently  small. 

Remark  1:  The  multidimensional  generalization  of  the 
computational  method  proposed  by  Kalmar-Nagy  and 
Huzmezan  involves  tracking  the  comer  points  of  the  boxes. 
The  comer  point  information  gives  the  edge  length.  The 
addition  rule  in  Eq  3  is  employed  for  each  linear  dimension 
which  in  turn  gives,  |Mfc+1|,  the  volume  of  the  box  in  (£+l)th 
step  around  xk+1  in  proportion  to  the  volume  of  the  box 
around  u.  Employing  Eq  7  in  the  iteration  eliminates  the 
need  to  keep  track  of  the  comer  points  since  the  expansion 
of  infinitesimal  volumes  under  the  map  can  be  computed  as 
the  determinant  of  Mk+l . 

IV.  Example 

We  now  illustrate  the  method  with  a  simple  example. 
Consider  the  feedback  system  with  G(x)  =  y/x .  The  pdf  of  u 
is  given  by  ....  Figure  2  shows  the  closed  loop  pdf  of  x  and 
the  iterates  of  the  pdf  converging  to  the  closed  loop  solution. 

V.  Parametric  uncertainty  in  the  loop  operator 

We  now  consider  the  case  when  additional  uncertainty  may 
be  present  in  the  loop  operator  G .  Suppose  that  G  depends 
on  a  parameter  0  whose  pdf  is  given  by  pg(0).  The  measure 
on  0  induces  a  measure  on  G,  pg-  As  in  Section  II,  we  can 
again  define  an  iteration  on  xk  and  Mk  with  the  difference 
being  that  G  is  drawn  randomly  at  each  iteration  step.  Let  Gk 
denote  the  sample  of  G  drawn  at  iteration  step  k  according  to 
pG.  Let  dGh  denote  the  derivative  of  Gk.  The  iteration  for  xh 
takes  the  form 

xk+l  =  u  +  Gkxk .  (8) 

This  system  belongs  to  the  class  of  iterated  random  func¬ 
tions  described  in  [4],  [5].  For  the  system  xk+1  =  f${xk) 


where  fg  is  a  random  map  selected  according  to  a  given 
probability  measure  p$9  under  the  assumption  that  fg  are  all 
Lipschitz  and  contractive  in  an  average  sense,  the  existence 
of  a  stationary  measure  and  the  convergence  are  described  in 
[4].  The  meaning  of  contraction  on  an  average  has  also  been 
quantified  in  terms  of  the  Lipschitz  constants  of  fg.  In  [5],  the 
contraction  on  average  assumption  is  relaxed  and  the  results 
are  extended  to  more  general  random  maps. 

VI.  Conclusion 

Iterative  methods  for  uncertainty  propagation  in  feedback 
systems,  together  with  graph  decomposition  methods  and 
block-by-block  propagation  can  lead  to  efficient  computational 
methods  for  complex  systems.  We  present  here  an  iterative 
scheme  to  compute  the  pdf  of  the  output  from  the  pdf  of  the 
input  for  a  nonlinear  feedback  loop.  The  convergence  of  the 
proposed  scheme  is  shown  under  the  assumption  that  the  loop 
operator  is  a  contraction.  Possible  extension  to  the  case  when 
there  is  additional  parametric  uncertainty  in  the  loop  operator 
is  described.  Future  work  will  focus  on  a  more  rigorous 
development  of  the  iterative  method  and  include  dynamical 
systems,  parametric  and  initial  condition  uncertainty  as  well 
as  integration  with  graph  decomposition  techniques  and  block- 
by-block  propagation  methods. 

Acknowledgments:  The  author  is  thankful  to  Tamas 
Kalmar-Nagy,  Mihai  Huzmezan,  Andrzej  Banszuk  at  UTRC 
and  Thordur  Runolfsson  at  University  of  Oklahoma  for  the 
insightful  discussions. 

References 

[1]  T.  Runolfsson,  I.  Mczic,  and  M.  Myers,  “Uncertainty  analysis  of  dynam¬ 
ical  systems,”  in  Proc.  of  the  SIAM  Conf.  on  Application  of  Dynamical 
Systems,  2003. 

[2]  S.  Varigonda,  T.  Kalmar-Nagy,  B.  LaBarrc,  and  I.  Mczic,  “Graph  de¬ 
composition  methods  for  complex,  interconnected  dynamical  systems,” 
in  Proc.  of  the  IEEE  CDC \  Bahamas,  2004,  (submitted). 

[3]  M.  Huzmezan  and  T.  Kalmar-Nagy,  “Propagation  of  uncertain  inputs 
through  networks  of  nonlinear  components,”  in  Proc.  of  the  IEEE  CDC, 
Bahamas,  December  2004,  (submitted). 

[4]  P.  Diaconis  and  D.  Freedman,  “Iterated  random  functions,”  SIAM  Review, 
vol.  41,  pp.  45-76,  1999. 

[5]  S.  F.  Jarncr  and  R.  Twcedie,  “Locally  contracting  iterated  functions  and 
stability  of  markov  chains,”  J.  Appl.  Prob.,  vol.  38,  pp.  494-507,  2001. 
[Online].  Available:  http://www.maths.lancs.ac.uk/jamer 

[6]  M.  Dcllnitz  and  O.  Junge,  Handbook  of  Dynamical  Systems  II:  Towards 
Applications .  World  Scientific,  2002,  ch.  Set  Oriented  Numerical 
Methods  for  Dynamical  Systems,  pp.  221-264.  [Online].  Available: 
http://math-www.upb.de/  agdellnitz/papcrs/handbook.pdf 

[7]  M.  Dcllnitz  and  R.  Preis,  Symbolic  and  Numerical  Scientific 
Computation ,  scr.  Lecture  Notes  in  Computer  Science,  2630. 
Springer,  2003,  ch.  Congestion  and  almost  invariant  sets  in 
dynamical  systems,  pp.  183-209.  [Online].  Available:  http://math- 
www.upb.de/  agdcllnitz/papcrs/congestion.ps.gz 

[8]  C.  Beck  et  al. ,  Thermodynamics  of  Chaotic  Systems .  Cambridge 
University  Press,  1995. 

[9]  T.-Y.  Li,  “Finite  approximation  of  the  Frobenius-Pcrron  operator,  a 
solution  to  Ulam’s  conjecture,”  J.  Approx.  Theory ,  vol.  17,  pp.  177-86, 
1976. 


Uncertainty  in  the  Dynamics  of  Conservative  Maps 


Oliver  Junge 
Institute  for  Mathematics 
University  of  Paderborn 
Warburger  Str.  100 
D-33098  Paderborn 
Germany 


Jerrold  E.  Marsden 
Control  and  Dynamical  Systems 
California  Institute  of  Technology,  MC  107-81 
Pasadena,  CA  91125 
USA 


Igor  Mezic 

Department  of  Mechanical 
and  Environmental  Engineering 
UCSB 

Santa  Barbara,  CA  93106 
USA 


Abstract — This  paper  studies  the  effect  of  uncertainty,  using 
random  perturbations,  on  area  preserving  maps  of  Wn  to  itself. 
We  focus  on  the  standard  map  and  a  discrete  Duffing  oscillator 
in  M2  as  specific  examples.  We  relate  the  level  of  uncertainty 
to  the  large  scale  features  in  the  dynamics  in  a  precise  way. 
We  also  study  the  effect  of  such  perturbations  on  bifurcations 
in  such  maps.  The  main  tools  used  for  these  investigations  is 
a  study  of  the  eigenfunction  and  eigenvalue  structure  of  the 
associated  Perron-Frobenius  operator  along  with  set  oriented 
methods  for  the  numerical  computations. 


I.  Introduction 


Uncertainty  is  obviously  of  central  importance  for  dy¬ 
namic  and  control  systems  both  from  a  theoretical  and 
a  practical  point  of  view.  While  this  topic  has  received 
much  attention  in  the  control  community,  there  is  still 
much  to  learn  about  the  effect  of  random  perturbations 
on  such  systems.  In  this  paper  we  focus  our  attention  on 
perturbations  of  conservative  system  as  a  first  step  towards 
studying  mechanical  systems  (including  molecular  systems) 
in  the  presence  of  uncertainty.  An  important  aspect  of  our 
approach  is  to  try  to  extract  the  key  dynamical  features  that 
survive  in  a  noisy  environment.  One  can  make  the  case  that 
such  features  are  the  most  important  ones  to  compute. 

The  specific  context  we  work  in  is  as  follows.  Let  /  : 
X  — >  X  be  an  area  preserving  map  of  a  compact  subset 
X  of  Rn  to  itself.  Let  B  be  the  Borel-cr-algebra  on  X 
and  p  the  standard  volume  measure  on  X.  We  model  a 
perturbed  version  of  the  map  /  by  a  stochastic  process 
which  maps  a  point  x  €  X  near  the  image  point  f(x) 
with  high  probability.  Formally,  for  a  given  perturbation 
size  6  >  0,  we  consider  a  stochastic  transition  function 
ps  :  X  x  B  — >  [0, 1]  corresponding  to  a  small  random 
perturbation  of  /  in  the  sense  of  [1],  i.e.,  we  require 

1)  x  »-»  p$(x,  B)  to  be  measurable  for  each  B  €  B; 

2)  ps(x,  •)  to  be  a  probability  measure  for  each  x  e  X; 

3)  for  all  continuous  functions  g  :  X  — >  R 


lim  sup  /  g(y)  ps{x,dy)  -  g{f(x)) 


=  0. 


(1) 


A  typical  small  random  perturbation  of  /  is  given  by  letting 
p${xi  •)  be  the  uniform  distribution  on  a  <5-ball  around  f(x ). 
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Throughout  the  paper,  we  focus  on  two  example  systems, 

1)  the  standard  map 

x\  \  _  /  x  +  ey  +  eVsin(27ra;), 

2/i  /  V  2/  +  e/i  sin(27rz) 

on  the  two-torus,  as  resulting  from  a  symplectic  time 
discretization  of  the  ordinary  differential  equation  for 
the  pendulum.  We  fix  the  step  size  e  =  1  and  view  p 
as  a  bifurcation  parameter; 

2)  and  the  Duffing  map , 

(  Xl  )  =  (  x  +  £y+£(px-x3),  \ 

\  2/i  /  V  V  +  f  (A4®  -  x3  +  pxi  -  xl)  J 

on  R2,  as  resulting  from  applying  the  Newmark 
scheme  to  the  Duffing  equations  (see,  e.g.,  [10]).  Here 
e  is  a  fixed  parameter  (corresponding  to  the  stepsize 
in  the  Newmark  scheme)  and  p  e  M  is  a  bifurcation 
parameter. 


A.  Stochastic  transition  functions 

A  measure  p  on  B  is  called  invariant  for  a  stochastic 
transition  function  p,  if 


/ 


p(x,A)  p{dx)  —  p(A) 


(4) 


for  all  A  e  B.  A  transition  function  p  :  X  x  B  — »  [0, 1]  is 
called  reversible  w.r.t  an  invariant  probability  measure  p,  if 

/  p(x,  B)  p(dx)  =  /  p(x,A)  p(dx).  (5) 
Ja  Jb 


The  value 


p(A  B)  '■=  y  lA  v{x,  B )  p(dx)  (6) 

is  the  probability  to  map  from  set  A  into  set  B  in  one  step. 
The  corresponding  probability  for  the  time  reversed  system 
is  given  by 

, M  m  p(B)p{B,  A) 

p{A’B)=  1%A)  •  (7) 


A  stochastic  transition  function  with  invariant  measure  p 
is  called  uniformly  ergodic ,  if  there  are  constants  q  <  1  and 
M  >  0,  such  that 


\\pn(x,-)-p\\<Mqn,  n  =  0,1,2,...,  (8) 


for  all  x  e  X. 

A  typical  small  random  perturbation  of  /  is  given  by 
letting  ps(x,  •)  be  the  uniform  distribution  on  a  (5-ball 
around  f(x).  This  is  an  example  of  an  absolutely  continuous 
transition  function ,  which  generally  is  of  the  form 

Ps(x,  B)  =  /  kg(f(x),y) 

JB 

where  k  :  X  x  X  — >  [0,  oo)  is  some  suitable  kernel.  In  this 
case,  the  transition  function  of  the  time  reversed  system  is 
explicitly  given  by 

P6(x,A)=  [  ks(f(y),x)dy(y).  (9) 

JA 

B.  Macroscopic  dynamics 

It  is  instructive  to  consider  the  situation  when  there  is 
no  noise,  for  an  area-preserving  mapping  T  (both  of  the 
maps  introduced  earlier  fall  into  this  category).  In  this  case, 
the  Perron-Frobenius  operator,  say  acting  on  functions  in 
L ?  is  given  by  Pf(x)  =  /oT"1  and  its  adjoint,  Koopman 
operator  by  Uf  =  foT.  These  are  both  examples  of  the  so- 
called  ’’transfer  operators”  ([9]).  The  spectral  properties  of 
these  can  be  related  to  phase-space  structure  of  conservative 
maps.  In  particular,  invariant  sets  of  T  can  be  obtained  as 
level  sets  of  functions  in  eigenspace  at  1  of  U  (see  e.g. 
[11]).  In  addition,  invariant  densities  of  T  are  functions  in 
eigenspace  at  1  of  P.  The  operators  P  and  U  are  unitary 
and  thus  their  spectrum  is  confined  to  the  unit  circle  in  the 
complex  plane.  Their  spectrum  can  be  analyzed  via  study 
of  a  simpler,  symmetrized  operator  ( U  +  P)/ 2.  Since  P 
is  normal,  any  eigenfunction  /  of  P  at  eigenvalue  A  is  an 
eigenfunction  of  U  at  eigenvalue  A,  the  complex  conjugate 
of  A.  In  that  case,  we  have 

1(1/  +  P)f  —  =  Re(A)/, 

and  thus  /  is  an  eigenvalue  of  the  symmetrized  operator 
(U  -j-  P)/ 2  associated  with  the  eigenvalue  Re(A).  In  the 
case  when  T  is  ergodic  (and  the  eigenspace  at  1  is  one¬ 
dimensional),  all  of  the  other  eigenspaces  are  at  most  two- 
dimensional.  However,  typical  area-preserving  maps  have  a 
more  complicated  structure  where  regular  regions  of  peri¬ 
odic  and  quasi-periodic  motion  coexist  with  chaotic  areas. 
An  example  of  this  is  shown  in  Figure  1,  for  the  standard 
map.  In  this  case,  the  eigenspace  at  1  can  be  obtained  using 
the  projection  operator  PTf  =  (1/n)  f  °  f.  An 
example  of  such  computation  is  given  in  figure  2. 

In  the  presence  of  noise  the  decomposition  of  the  state 
space  of  an  area-preserving  map  into  invariant  sets  with 
“regular”  resp.  “chaotic”  dynamics,  as  described  above,  (cf. 
Figure  1)  is  no  longer  possible. 

Due  to  the  noise,  the  evolution  may  become  transitive  on 
the  state  space.  However,  for  certain  types  of  perturbations 
and  if  the  perturbation  is  sufficiently  small,  certain  invariant 
sets  of  the  unperturbed  map  can  still  be  recovered  as  almost 
invariant  subsets.  In  fact,  this  is  the  more  likely  thing  one 


Fig.  1.  Standard  map:  500  iterates  of  500  different  initial  conditions, 
chosen  at  random  according  to  a  uniform  distribution,  y,  =  0.3. 


X 


Fig.  2.  Standard  map:  Invariant  sets  obtained  from  Prf,  for  /  — 
cos(27ry),^  —  0.3. 


will  see  in  a  more  realistic  noisy  and  uncertain  model  of  a 
system.  Intuitively,  a  set  A  C  X  is  almost  invariant ,  if  the 
invariance  ratio  of  A , 

P6(A,A)9  (10) 

is  close  to  1,  i.e.  if  the  probability  to  stay  within  the 
set  A  under  the  evolution  is  high.  We  refer  to  [2]  for 
a  more  formal  definition.  For  the  perturbed  system,  it  is 
therefore  natural  to  ask  for  a  decomposition  of  the  state 
space  into  almost  invariant  sets  when  one  is  interested  in  a 
more  macroscopic  and  robust  description  of  the  dynamical 
behavior.  Thus,  we  aim  at  finding  sets  A  for  which  (10)  is 
as  large  as  possible.  This  question  has  actually  been  treated 
in  [3],  [4]  in  detail.  It  is  a  trivial  (but  important)  observation, 
that 

Ps(A,A )  =ps(A,A)  =  ±(ps{A,A)  +MAA)),  (11) 


i.e.  when  looking  for  sets  that  maximize  ps{A ,  A),  we  might 
as  well  consider  the  transition  function 

r6(x,A)  =  ^(ps(x,A)  +  p6(x,A)),  (12) 

which  is  reversible. 

C.  Transfer  operators 

As  mentioned  before,  spectral  properties  of  transfer  oper¬ 
ators  can  be  used  to  explore  transport  in  phase  space.  One 
possibility  to  identify  sets  with  a  large  almost  invariance  ra¬ 
tio  is  via  a  spectral  analysis  of  a  transfer  operator  associated 
to  the  evolution  law  under  consideration,  see  e.g.  [2].  For 
a  given  (“pointwise”)  evolution  law,  an  associated  transfer 
operator  describes  the  evolution  of  probability  measures  or 
densities  on  the  state  space.  For  a  deterministic  map  /,  the 
Perron-Frobenius  operator  P  :  M  — >  M, 

Pp(A)  =  p(/_1(A)),  AeB ,  (13) 

where,  instead  of  L2(M)  here  P  is  defined  on  M,  the 
space  of  probability  measures  on  ( X ,  B).  The  corresponding 
transfer  operator  of  a  stochastic  transition  function  p$  is 
given  by  Ps  :  M  — >  M, 

Psn(A)  =  J ps(x,A)  dfj,(x).  (14) 

In  the  case  that  ps  is  absolutely  continuous,  we  can 
equivalently  consider  Ps  as  an  operator  on  L2(X)i  and 
Ps  :  L2  — ►  L 2  is  given  by 

( Psh)(y )  =  J  kg(f(x),y )  h(x)  dm(x). 

Note  that  for  the  standard  inner  product  on  L2,  the  adjoint 
of  Ps  is  given  by 


n  is  P(m  =  l  !%  f(0-  w  -  0 df.  If  we  set  f{6)  = 

exp{i2icn6),  we  obtain 

1  rsI 2 

P/(0)  =  T  /  exp(i27rn(0  -  a;  -  f))<*C 

d  J-(5/2 

=  i  exp(— i27muj)  exp(i27cn9) 
o 

,6/2 

/  exp(— z27rn£)(i£ 

J-8/2 

_  sin (mr5)  rnu)  exp(i27m0) 


Therefore,  exp(i2nn0)  are  eigenfunctions  of  P%  associated 
with  eigenvalue  An  =  exp(— i27muj).  For  fixed  n, 

if  J  — >  0,  An  — *  1.  In  figure  3  we  show  eigenvalues  of  P^ 
for  u)  =  7t/320,  6  =  0.01.  For  large  n  the  eigenvalues  tend 
to  zero.  The  eigenvalue  with  the  largest  modulus  smaller 
than  1  is  Ai  =  exp(— i27rtu),  with  the  associated 

eigenfunction  exp(i27T0).  The  adjoint  operator  of  P$  is  the 


Fig.  3.  Eigenvalues  of  for  u)  =  7r/320,  <5  =  0.01  shown  as  full  solid 
circles  inside  the  unit  circle 


(pl9)(x)  =  J ks(f(x),y)  g(y)  dm(y). 

We  first  study  an  example  of  the  Perron-Frobenius  op¬ 
erator  for  a  stochastic  perturbation  of  a  conservative  map 
on  a  circle;  specifically,  we  discuss  the  spectral  properties 
of  the  Perron-Frobenius  operator  associated  with  rotations 
on  the  circle.  In  this  example,  spectral  properties  can 
be  analytically  computed  and  as  a  result  are  simple  to 
understand.  Let  T  :  Sl  — >  S1;#'  =  0  +  u,  be  a 
rotation  on  the  circle,  where  a;  is  a  constant.  Note  that 
this  map  preserves  the  Haar  measure  on  the  circle.  The 
Perron-Frobenius  operator  associated  with  it  is  given  by 
Pf(0)  =  The  eigenfunctions  of  P  associated  with 

eigenvalue  exp(-227 mu;),  n  £  N  are  easily  computed  to  be 
exp(i27rn$).  These  eigenfunctions  span  L2{Sl).  Consider 
the  following  stochastic  perturbation  of  this  conservative 
system:  :  6'  =  6  +  u)  +  £,  where  f  is  a  random 

variable  uniformly  distributed  on  [—5/ 2,tf/2],  where  (5  is 
a  constant.  The  Perron-Frobenius  operator  associated  with 


Koopman  operator  |  f-s/2  /(^  +  ^ 

this  case  the  Perron-Frobenius  and  Koopman  operators  are 
normal  and  thus  the  symmetrized  operator  (P^  +  Uf)/2 
has  eigenvalues  (sm(n7r8)(n7t6)  cos(27rnca))  and  the  same 
eigenfunctions  as  For  small  u  the  second  largest  eigen¬ 
value  is  obtained  for  n  =  1.  The  associated  eigenfunctions 
are  given  by  c\  sin(27r$)  +  c<i  cos(27T0).  Clearly  any  two 
arcs  that  split  the  circle  in  two  are  almost  invariant  sets. 

For  a  reversible  and  uniformly  ergodic  stochastic  transi¬ 
tion  function  p,  the  essential  spectral  radius  of  the  corre¬ 
sponding  transfer  operator  P  :  L2  — >  L2  is  bounded  away 
from  one,  i.e.  cr(P)  C  r]  U  {1}  with  r  <  1;  see  [6]. 
By  combining  Theorem  3.1  and  Corollary  4.33  of  [6],  we 
obtain: 


Proposition  1  Let  the  stochastic  transition  function  p  be 
reversible  and  uniformly  ergodic.  Assume  that 

a(P)c  MU{A2)1}, 


(16) 


with  r  <  X2  and  X2  a  simple  eigenvalue  of  P,  Then  for  any 
set  A  £  A 

l-\2<p{A,Ac)+p{Ac,A)<l-K\2,  (17) 

where  0  <  k  <  1. 

In  fact,  one  has  (see  [6]) 


where  v2  is  the  eigenvector  of  P  corresponding  to  the 
eigenvalue  A,  i.e.  k  is  measuring  the  “L2-deviation”  of  v2 

from  ffff-XA  -  So,  roughly  speaking,  the 

maximal  almost  invariance  ratio  of  any  set  A  is  given  by 
(A2  +  l)/2. 

Also,  the  form  of  k  suggests  a  strategy  for  the  identi¬ 
fication  of  A,  given  v2:  the  quantity  will  be  maximized, 
if  the  set  A  is  defined  to  be  the  subset  of  X ,  on  which 
v2  is  positive.  This  strategy  of  identifying  almost  invariant 
subsets  of  phase  space  has  successfully  been  applied  to  the 
identification  of  conformations  of  molecules,  see  [5],  [12]. 

In  Figure  4  we  show  six  eigenvectors  of  the  reversibilized 
discretized  transfer  operator  (of  the  unperturbed  system). 
Clearly  the  almost  invariant  sets  defined  by  the  positive 
resp.  negative  components  of  these  modes  are  in  very  good 
agreement  with  invariant  sets  of  the  unperturbed  system  (c.f. 
Figure  1)  for  the  eigenvectors  v2i . . .  ,v6.  The  eigenvector 
V7  decomposes  the  “chaotic  sea”  itself. 

D.  Discretization 


(a)  v2 


(b)  *>3 


I 


(c)  Vi 


1 


_ J 


(C)  V6 


(f)  V7 


Numerically,  we  need  to  work  with  a  finite  dimensional 
subspace  Ld  of  L 2  and  a  corresponding  approximation  of 
P  on  Ld-  A  common  ansatz  is  to  let  Ld  be  the  subspace  of 
step-functions  that  are  piecewise  constant  on  the  elements 
of  a  partition  of  X.  More  precisely,  let  V  =  {B\, . . . ,  Bd } 
be  a  finite  subset  of  B  such  that  p(Bi  n  Bj)  =  0  for 
Bi,Bj  €  V,Bi  ^  Bj  and  ULi  Bi  =  x-  Let  L*  = 
span{xu1,  ■  • .  ,XBd}  and  let  Qd  :  L2  -+  Ld  be  the  orthog- 
onal  projection  onto  Ld-  The  discretized  transfer  operator 
Pd  is  defined  as 

Pd  =  QdPs •  (18) 


It  is  easy  to  see  (see  e.g.  [2])  that  in  this  case  the  discretized 
operator  can  be  represented  by  a  stochastic  matrix  with 
entries 


™(/-Wngj) 

Pl]  m{Bj) 


i,j  = 


(19) 


For  bounded  perturbations  of  the  underlying  map  /,  i.e.  if 
p(Xj-)  has  sufficiently  small  support,  then  this  transition 
matrix  is  sparse  and  even  if  the  number  d  of  partition 
elements  is  quite  large,  a  couple  of  eigenvalues  and  cor¬ 
responding  eigenvectors  can  efficiently  be  computed  by 
Amoldi  methods  as  e.g.  implemented  in  ARPACK  [8] 
(which  is  available  as  the  eigs  command  in  Matlab). 


Fig.  4.  Standard  map:  eigenvectors  V2,---,V7  to  the  six  largest 
eigenvalues  (except  1)  A2  <  •  •  •  <  A7  of  the  reversibilized  discretized 
transfer  operator  of  the  unperturbed  system  on  a  partition  of  212  boxes, 
H  =  0.3. 


As  an  example,  again  we  have  a  look  at  the  standard  map. 
We  consider  the  reversibilized  discretized  transfer  operator, 
which  in  this  case  is  given  by  the  matrix  Rd  =  (Pd+Pj)/ 2 
and  the  decomposition  of  the  state  space  given  by  the  sign 
structure  of  the  eigenvector  corresponding  to  the  second 
largest  eigenvalue.  In  Figure  5  we  show  this  eigenvector, 
computed  on  four  different  partitions  of  the  phase  space. 
Note  that  the  overall  shape  of  the  decomposition  is  already 
quite  well  resolved  on  a  partition  with  only  256  boxes. 

II.  Spectral  Properties  of  the  Transfer 
Operator 

A.  The  dependence  of  the  spectrum  on  the  perturbation  size 

Let  us  focus  on  a  particular  perturbation  of  the  given  map 
/:  We  consider 

kg  {f(x),y)  =  — Xbs  (/(*))  (y)  (20) 

In  [7]  we  prove: 


(c)  2 12  boxes 


(d)  216  boxes 


Fig.  5.  Standard  map:  eigenvector  to  the  second  largest  eigenvalue  of  the 
reversibilized  discretized  transfer  operator  of  the  unperturbed  system  for 
different  partitions,  ^  =  0.3. 


Proposition  2  If  A  is  an  invariant  set  of  the  unperturbed 
system,  then  the  probability  p$(A,  Ac)  to  map  from  A  into 
its  complement  Ac  under  the  perturbed  system  can  be 
estimated  as 

Ps(A,  Ac )  <  lCTff  - S  +  0(S2)  (21) 

as  6  — *  0. 

Combining  this  estimate  with  Proposition  1,  we  obtain 


Theorem  1  Let  A  be  an  invariant  set  of  the  unperturbed 
map  f.  Assume  that  the  reversibilized  small  random  per¬ 
turbation  rs  —  \{p5  +ps)  of  f  is  uniformly  ergodic .  Then 


1-A2  < 


length(dA) 

MW) 


5  +  0(52), 


(22) 


Proof:  Note  that  p{A)-\~p{Ac)  =  1  and  length (5 A)  = 
length(dAc).  ■ 

Note  that  the  constant  in  c 2  relates  the  length  of  the 
boundary  of  A  to  its  volume.  In  particular  this  means  that 
a  larger  eigenvalue  corresponds  to  an  almost  invariant  set 
of  larger  volume.  So  in  this  sense  larger  eigenvalues  detect 
“more  important”  almost  invariant  sets. 

In  Figure  II- A  we  show  how  the  four  largest  real  eigen¬ 
values  of  the  discretized  transfer  operator  (14  subdivisions, 
214  boxes)  for  the  standard  map  depend  on  the  perturbation 
size  5.  For  the  computation  of  the  transition  matrix  we  used 
16  test  points  in  a  regular  grid  in  each  box  of  the  covering 
and  256  points  on  a  regular  grid  in  each  ^-neighborhood  of 
the  corresponding  image  points  in  order  to  sample  ps(x,  •). 


Note  that  this  numerical  result  is  in  very  good  agreement 
with  Theorem  1. 


Standard  map,  eigenvalues  vs  perturbation  size,  depth  14. 16f256  Grld/InnerGrid 


Fig.  6.  Standard  map:  the  four  largest  real  eigenvalues  (except  1)  in 
dependence  of  the  perturbation  size. 


B.  Spectrum  for  a  fixed  value  of  the  bifurcation  parameter 

Suppose  that  we  are  given  a  symmetry  transformation 
n  :  X  — >  X,  such  that  n2  =  id  and  the  map  /  satisfies  the 
equivariance  condition 

no  f  —  /_1  o  n.  (23) 

Define  the  action  of  n  :  L2  —>  L2  by 

nh  =  ho  n. 


Proposition  3  Suppose  that  the  kernel  k  satisfies 


h{f(x),y)  =  ks(Ko  f  '(y),K{ x)), 

(24) 

then 

nPs  =  Pi  n. 

(25) 

Proof:  We  compute 

pj Kh(x)  =  J  ks(f(x),y)  h  o  n(y)  dy(y) 

(26) 

=  J  ks(K0f  1(y),n(x))  ho  K{y)  dfi(y)  (27) 

=  J ks(K(y'),K(x))  ho  ko  f(y')  dy(y'),  (28) 

where  the  latter  equality  follows  from  the  change  of  vari¬ 
ables  y*  =  f~1(y)  and  the  fact  that  /  is  area  preserving. 
Using  the  equivariance  condition  (23)  we  continue 


=  [  k5(n(y'),K(x))  hof  1  o  K(y')  dy{y') 

u 

(29) 

=  J  ks(y,  k(x))  h  o  f~\y)  dy(y) 

(30) 

=  J  ks(f{y),K(x))  h{y)  dy(y) 

(31) 

=  P$h  o  n(x)  =  nPsh(x ), 

(32) 

where,  again,  we  performed  two  changes  of  variables, 
exploiting  the  area  preservation  property  of  n.  ■ 

Corollary  1  Let  h  e  L2  be  an  eigenfunction  of  Ps  to  the 
eigenvalue  A,  then  nh  is  an  eigenfunction  of  Pj  to  the 
eigenvalue  A. 

Numerical  support  for  this  Corollary  is  shown  in  Figure  7. 
We  computed  the  transition  matrix  on  a  covering  of  212 
boxes,  using  100  points  on  a  regular  grid — without  any 
perturbation  (i.e.,  S  =  0),  and  for  e  =  0.1  and  p  =  —0.2. 
The  figure  shows  ei7Tnh  in  the  left  column,  where  h  is 
the  eigenvector  of  Ps  and  the  eigenvector  of  Pj  in  the 
right  column,  both  for  the  same  eigenvalue  A  =  0.9902  4- 
i  0.0647. 
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Fig.  7.  Duffing  map:  Comparison  of  el7r  nh  (left  column),  where  h  is 
the  eigenvector  of  Ps,  and  the  corresponding  eigenvector  of  Pj  (right 
column). 


C.  Spectrum  in  dependence  of  the  bifurcation  parameter 

As  a  bifurcation  parameter  of  the  map  under  consider¬ 
ation  is  varied,  it  is  natural  to  expect  the  spectrum  of  the 
(discretized)  transfer  operator  to  change.  In  particular,  we 
expect  the  corresponding  eigenmodes  to  reflect  changes  in 
the  global  dynamical  behavior  of  the  system.  As  an  example 
situation  we  consider  the  dynamics  of  the  Duffing  map  as 
the  parameter  p  is  varied.  It  is  well  known  that  at  p  —  0 
the  origin  undergoes  a  pitchfork  bifurcation,  leading  to  two 
stable  equilibria  for  p  >  0.  In  Figure  8  we  show  how  part 
of  the  spectrum  of  the  reversibilized  discretized  transfer 
operator  depends  on  p  £  [-0.35,0.7]. 


One  particular  eigenfunction  (associated  to  the  eigenvalue 
represented  by  the  red  curve  in  Figure  8)  appears  to  become 
the  dominant  one  as  p  changes  from  -0.35  to  0.7.  For 
positive  p,  this  positive  resp.  negative  components  of  this 
eigenfunction  are  associated  to  neighborhoods  of  the  two 
equilibria  which  are  created  in  the  pitchfork  bifurcation, 
bounded  by  the  homoclinic  orbits  of  the  origin.  For  negative 
p  however,  this  geometric  situation  is  not  present,  but  still 
the  sign  structure  of  this  eigenmode  defines  a  qualitatively 
similar  decomposition  of  the  state  space  into  two  almost 
invariant  sets. 


III.  Conclusions. 

In  this  paper  we  have  shown  that  one  can  effectively 
use  the  theory  and  computation  associated  with  the  Perron- 
Frobenius  operator  to  study  the  effect  of  uncertainty  on 
area  preserving  maps  and  illustrated  the  methods  with  the 
standard  map  and  a  discrete  Duffing  oscillator.  The  level  of 
uncertainty  is  related  in  a  quantitative  way  to  the  large  scale 
features,  often  the  ones  that  are  the  most  important  to  first 
compute.  It  was  also  shown  the  way  in  which  uncertainty 
affects  the  bifurcations  of  such  maps. 
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Abstract. 

This  report  describes  the  form  of  the  dynamical  systems  involved  in  computational  chemistry  and 
molecular  modelling  and  points  out  the  avenues  where  methods  of  uncertainty  analysis  and  uncertainty 
propagation  may  be  usefully  applied.  Specific  possibilities  include  the  quantification  of  uncertainty  in  the 
semi-empirical  potential  energy  functions  in  classical  molecular  dynamics  and  the  assessment  of  uncertainty 
in  quantum  molecular  dynamics  (Car-Parrinello  molecular  dynamics)  due  to  the  incomplete  knowledge  of 
the  exchange  correlation  functional.  These  can  have  bearings  on  atomistic  simulation  in  general.  In  relation 
to  this,  two  example  dynamical  system  definitions  are  given  which  may  be  starting  points  in  the 
investigation  of  the  Dellnitz-Preis  and  polynomial  chaos  methods  in  this  area. 

1.  Introduction. 

Experimental  values  of  physical  quantities  of  a  many-particle  system  can  be  found  as  an  ensemble 
average.  Experimental  systems  are  so  large  that  it  is  impossible  to  determine  this  ensemble  average  by  sum¬ 
ming  over  all  accessible  states  on  a  computer.  There  exist  essentially  two  methods  for  determining  these 
physical  quantities  as  statistical  averages  over  a  restricted  set  of  states:  the  molecular  dynamics  (MD)  and 
Monte  Carlo  (MC)  methods.  Suppose  we  have  a  random  sample  of,  say,  107  configurations  of  the  system 
which  are  all  compatible  with  the  values  of  the  system  parameters.  For  such  a  large  number  we  expect  aver¬ 
ages  of  physical  quantities  over  the  sample  to  be  rather  close  to  the  ensemble  average.  It  is  unfortunately 
impossible  to  generate  such  a  random  sample;  however,  we  can  generate  a  sample  consisting  of  a  large  num¬ 
ber  of  configurations  which  are  determined  successively  from  each  other  and  are  hence  correlated.  This  is 
done  in  molecular  dynamics  and  Monte  Carlo  methods.  Our  focus  in  this  report  is  on  molecular  dynamics. 

Molecular  dynamics  is  a  widely  used  method  for  studying  classical  many-particle  systems.  It  con¬ 
sists  essentially  of  integrating  the  equations  of  motion  of  the  system  numerically.  It  can  therefore  be  viewed 
as  a  simulation  of  the  chemical  or  molecular  system  as  it  develops  over  a  period  of  time.  The  system  moves 
therefore  in  phase  space  along  its  physical  trajectory  as  determined  by  the  equations  of  motion.  The  great 
advantage  of  MD  is  that  it  not  only  provides  a  way  to  evaluate  expectation  values  of  static  physical  quanti¬ 
ties;  dynamical  phenomena  such  as  transport  of  heat  or  charge,  or  relaxation  of  systems  far  from  equilibrium 
can  also  be  studied. 

The  remainder  of  this  report  is  organised  as  follows.  Section  2  describes  classical  molecular 
dynamics  and  points  out  avenues  for  uncertainty  analysis  in  such  calculations.  Section  3  describes  quantum 
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molecular  dynamics  in  the  context  of  atoms  in  molecules.  The  verlet  algorithm  for  classical  and  quantum 
molecular  dynamics  simulations,  a  symplectic  integrator,  is  described  for  the  situation  with  holonomic  con¬ 
straints.  The  importance  of  uncertainty  analysis  due  to  the  inexact  knowledge  of  the  exchange  correlation 
functional  is  emphasised  a  the  most  important  problem  to  be  tackled  either  by  polynomial  chaos  or  any  other 
method  of  uncertainty  propagation.  Section  4  describes  molecular  systems,  in  particular  the  nitrogen  mole¬ 
cule  and  the  CS2  trimer,  which  can  be  used  as  model  cases  to  study  different  methods  of  uncertainty  analysis 
-  Dellnitz-Preis  graph  theoretic  methods  or  stochastic  galerkin  methods. 


2.  Classical  Molecular  Dynamics. 

Consider  a  collection  of  N  classical  particles.  The  particles  interact  with  each  other  via  an  interac¬ 
tion  potential  which  is  the  negative  gradient  of  a  force-field.  If  the  internal  force  acting  on  particle  i  be 
denoted  by  F^R)  where  R  comprises  the  position  coordinates  rx  of  all  particles.  In  the  absence  of  gravita¬ 
tional  forces  and  other  forces  due  to  boundaries,  the  equations  of  molecular  dynamics  are 


n(t) 


Fj(R) 


(1) 


in  which  m,-  is  the  mass  of  particle  i.  The  solutions  of  the  equations  of  motion  describe  the  time  evolution  of 
a  real  system  although  the  description  of  equation  (1)  is  approximate  for  the  following  reasons: 

(a.)  It  is  a  classical  description  as  per  Newton’s  laws  of  motion.  The  importance  of  quantum  effects  depends 
strongly  on  the  particular  type  of  system  considered  and  on  the  physical  parameters  (temperature,  density, 
etc.).  In  Section  3,  quantum  molecular  dynamics  is  described  in  the  context  of  density  functional  theory. 

(b.)  The  forces  between  the  particles  are  not  known  exactly;  quantum  mechanical  calculations  from  which 
they  can  be  determined  are  subject  to  systematic  errors  as  a  result  of  the  neglect  of  correlation  effects,  as  we 
have  seen  in  previous  chapters.Usually  these  forces  are  given  in  a  parametrised  form,  and  the  parameters  are 
determined  either  by  ah  initio  calculations  or  by  fitting  the  results  of  simulations  to  experimental  data. 

In  a  molecular  dynamics  simulation,  one  initialises  a  system,  starts  the  simulation,  and  lets  the  sys¬ 
tem  reach  equilibrium.  The  number  of  particles  and  the  form  of  the  interaction  are  specified.  The  particles 
are  assigned  positions  and  momenta.  If  a  Lennard-Jones  potential,  which  is  discussed  below,  is  used,  the 
positions  are  usually  chosen  as  the  sites  of  a  Bravais  face  centered  cubic  (fee)  lattice  [1],  which  is  the  ground 
state  configuration  of  the  Lennard-Jones  system.  The  velocities  are  drawn  from  a  Maxwell  distribution  at  the 
specified  temperature.  The  numerical  solution  of  molecular  dynamics  is  achieved  by  symplectic  integrators 
such  as  the  Verlet  algorithm  [2].  In  Section  3.1  a  description  of  the  Verlet  method  in  the  presence  of  holo¬ 
nomic  constraints  is  given. 
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The  development  of  force  fields  is  based  on  a  judicious  combination  of  accurate  electronic  struc¬ 
ture  calculations  and  some  empiricism.  The  empiricism  in  the  force  fields  results  in  some  uncertainty  whose 
propagation  in  the  molecular  dynamics  simulation  needs  to  be  quantified.  Over  the  years,  a  number  of  poten¬ 
tials  such  as  bond-order  potentials,  embedded  atom  potentials,  and  simple  pair  potentials  have  been  devel¬ 
oped;  as,  e.g.,  see  equations  (5)  -  (10)  in  reference  [3];  a  full  description  of  empirical  potentials  may  be 
found  in  Chapter  4  of  reference  [4].  A  famous  force  field  is  the  Stillinger  Weber  potential  used  in  the  molec¬ 
ular  dynamics  calculations  of  Silicon  and  other  semiconducting  systems.  The  parameters  in  these  force 
fields  are  determined  by  fitting  to  experimental  structure,  phonon  spectra,  and  mechanical  properties.  Uncer¬ 
tainty  in  these  fitted  parameters  can  have  an  impact  on  the  computed  properties  obtained  from  simulating 
equation  (1).  The  use  of  polynomial  chaos  in  order  to  assess  the  propagation  of  uncertainty  in  molecular 
dynamics  with  respect  to  the  uncertain  parameters  in  the  interatomic  potential  of  force-fields  is  an  as  yet  un¬ 
investigated  research  problem. 

An  example  of  an  interatomic  force-field  is  the  negative  gradient  of  a  van  der  Waals  potential  func¬ 
tion,  the  Lennard- Jones  12-6  function.  The  Lennard-Jones  12-6  function  takes  the  following  form 


(2) 


where  the  two  parameters  are  the  collision  diameter  a  and  the  well-depth  e.  The  effect  of  uncertainty  in  a 
and  e  on  molecular  dynamics  is  worthy  of  investigation  by  polynomial  chaos  methods.  Another  challenging 
problem  is  the  application  of  polynomial  chaos  methods  to  molecular  dynamics  in  the  context  of  the  Still¬ 
inger- W  eber  potential;  see  equations  (4.121)  -  (4.123)  in  reference  [4]. 


3.  Quantum  Molecular  Dynamics:  Atoms  in  Molecules. 

Given  a  diatomic  molecule  AB  and  its  ground  state  electron  density  p AB(r)  obtained  by  the  stan¬ 
dard  Hohenberg-Kohn-Sham  density  functional  theory  (DFT)  [5],  we  seek  to  solve  the  optimisation  problem 

minimise  EA(pA;RA)  +  EB(pB;RB) ,  (3) 

subject  to 

P  abM  =  P,^)  +  Pb(')  Vr,  (4) 

in  order  to  determine  pA(r)  and  p B(r) ;  RA  and  R&  are  held  fixed,  and  are  the  respective  nuclear  coordinates 
of  the  atoms  A  and  B  for  a  given  intemuclear  separation.  EA  and  EB  are  the  total  energy  functionals  of  atoms 

A  and  B.  In  terms  of  the  effective  single-particle  orthonormal  orbitals  {\|/f  (r)>  and  {\}/f  (r)}  , 
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and 


p^')  =  X«ik*(r)|2 

k 


(5) 


k 


(6) 


where  and  nf  are  the  occupation  numbers  for  the  one-particle  wave-functions  \|/*(/*)  and  \j/£(r)  respec¬ 
tively.  The  Kohn-Sham  total  energy  functional  for  each  isolated  atom  is  given  by 


EAUvf(r)},*A)*EA{pA,RA)  =  I«*J(V^»-))*[-^V2]^(r)rfi--e2z4^Vr+  (7) 


2ffP^(»,)Pj4(0 


2JJ 


|r-r'| 


tfrA-'  +  S^p^r)) 


and 


**((¥?(»■)}.  *B)-S*(P»*a)  =  S"*J('l,f('-))*[-8^V2]vl/f(r)^-e2ZflJ^V,+ 


(8) 


2  f  f  P^(r)P  B^r  )j  j,,  P  /  ^ 

e  JJ  TT 7i — +  £*c(Ps(r)) 


|r-  r'| 


where  e  is  the  electronic  charge,  h  is  Planck’s  constant,  m  the  electronic  mass,  and  ZA  and  ZB  are  the  respec¬ 
tive  atomic  numbers  oft4  and  B.  Exc  is  the  exchange  correlation  functional.  The  spin  orbitals  must  satisfy  the 
orthonormality  constraints 


<ViW>  =  hi  (9) 

and 

<vflvf>  =  8w.  (10) 

The  solution  of  the  optimisation  problem  given  by  equations  (2)  and  (4)  may  be  achieved  by  quan¬ 
tum  molecular  dynamics  via  a  slight  modification  of  the  dynamical  simulated  annealing  method  of  Car  and 
Parrinello  [6].  This  requires  the  construction  of  a  dynamical  Lagrangian  which  includes  the  electronic  wave 
functions  and  their  time  derivatives  as  the  variables  with  respect  to  a  fixed  internuclear  separation.  This  leads 
to  a  classical  mechanics  problem  with  the  sum  of  the  energy  functionals  of  equations  (7)  and  (8)  acting  as  a 
potential.  If  a  friction  term  is  then  added  to  the  equations  of  motion  of  this  classical  system,  the  degrees  of 
freedom  will  come  to  rest  after  some  time,  with  values  corresponding  to  the  minimum  of  the  classical  poten¬ 
tial,  which  is  the  energy  of  the  quantum  system  at  the  equilibrium  configuration  of  the  nuclei.  It  is  also  pos- 
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sible  to  set  the  frictional  force  to  zero  in  order  to  simulate  the  system  at  nonzero  temperature.  The 
Lagrangian  is  given  by 


(r)},  {yf (r)}^,  rb)  =  y-X^t)2  +  yX(vf)2 -  (1 1) 

k  k 

EA{{^{r)}-RA)  -  £s({y?  (I-)}*,)  +  £2>wKV*lW>  -  8W]  + 

k  l 

XXAW[<'hW?>  -  8tf]  +  f  EAB(r)[pA(r)  +  p B(r)  -  p AB{r)]dr 
k  l 

where  ) iA  and  \xB  are  small  masses  and  MA  and  MB  are  the  respective  masses  of  the  nuclei  of  A  and  B  respec¬ 
tively.  ,  hBkl ,  and  V AB(r)  are  the  Lagrange  multipliers  which  enforce  orthonormality  and  the  atoms  in 

molecule  constraint  of  equation  (4).  There  is  no  term  in  the  Lagrangian  of  equation  (11)  that  pertains  to  Cou- 
lombic  repulsion  between  the  nuclei  of  A  and  B  because  it  is  a  constant  at  fixed  internuclear  separation.  The 
details  of  the  kinetic  energy  of  the  electrons  do  not  matter;  the  fictitious  masses  ]XA  and  paused  in  the  kinetic 
energy  term  can  be  set  to  unity.  If  friction  is  included  into  the  equations  of  motion,  the  particular  values  of 
the  electronic  and  nuclear  masses  do  not  matter  as  the  kinetic  energy  is  zero  at  the  minimum  of  the  potential 
of  the  Lagrangian  of  equation  (11).  The  Euler-Lagrange  equations  corresponding  to  the  first  variations  of  the 
functional  of  equation  (1 1)  are  given  by 


VAVk  = 


:EAttV?(r)},  Ra)  +  £Ah(OM ,f(r)  +  J VAB(r) 


dr 


(12) 


and 


wt 


0\|/‘ 


;EBa  Vf(r)},  Rb)  +  £A%mf  (r)  +  J  V AB(r) 


AM') 


dr 


(13) 


The  time  integration  of  the  quantum  molecular  dynamics  equations  (12)  and  (13)  can  be  undertaken  by  the 
use  of  Verlet’s  method. 

3.1  Verlet’s  Method  for  Constrained  Dynamics. 

Verlet’s  method  is  widely  used  for  the  simulation  of  many-particle  dynamics  in  the  computer  simu¬ 
lation  of  liquids  [2].  Suppose  there  are  N  particles  with  R  denoting  the  position  coordinates  rt  of  all  the  parti¬ 
cles.  Consider  the  dynamics 
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j2 

-^-r.(0  =  Ff(*)  (14) 

dt 

where  is  a  force-field.  Verlet’s  method  in  the  absence  of  any  constraints  on  the  particles  is  the  recurrence 
relation 


r{t  +  h)  =  2rt(t)  -  rt(t-  h)  +  h2Fj[r(t)]  (15) 

which  is  accurate  to  0(h4)  and  wherein  r-(/)  is  the  position  of  the  ith  particle  at  time  t  =  nh  where  h  is  the 
time-step  and  n  is  an  integer.  Suppose  M  constraints,  ck(R)  =  0,  k  =  1, M ,  are  imposed  upon  the  parti¬ 
cles.  The  constraint  force  acting  on  the  particle  is 

M 

X  h^rP(R)  (16) 

1 

where  the  {^}  are  the  Lagrange  multipliers  to  be  determined.  At  time  t  =  nh  we  have  at  our  disposal  the 
positions  at  times  t  and  t-h  .  First,  one  calculates  new  positions  r^t  +  h)  ignoring  the  constraints  according 
to 


r,(f  +  A)  =  2r<(/)-ri(r-A)  +  A2/?l[r(/)]  (17) 

and  the  corrected  new  positions  are  given  by 

M 

r^t  +  h)  =  rfif  +  h)-  X  XkVrp(R)  (18) 

k=  1 

where  the  are  found  by  an  iterative  procedure.  The  iterations  are  numbered  by  an  index  /.  In  each  itera¬ 
tion,  a  loop  over  the  constraints  k  is  carried  out,  and  in  each  step  of  this  loop,  the  Lagrange  multiplier  X^  and 
all  the  particle  positions  are  updated.  The  positions  are  updated  according  to 

r"ew  =  rf-h2x"\rpk(R(t)).  (19) 

The  parameter  X*  *  is  found  from  a  first  order  expansion  of  ok(R(t))  and  requiring  that  the  latter  vanishes: 

ct(Rne"')  =  ak(R°ld)  -  '£Vrpk(Rold)Vrpk(R(t))  =  0 .  (20) 


This  leads  to 


Each  step  will  therefore  shift  the  positions  more  closely  to  the  point  where  they  all  satisfy  the  constraint.  The 
iterations  are  stopped  when  all  the  constraints  are  smaller  in  absolute  value  than  some  small  positive  number. 

3.2  Investigation  of  Uncertainty  in  the  Exchange  Correlation. 

Using  polynomial  chaos  to  study  the  propagation  of  uncertainty  in  Car-Parrinello  calculations  is  a 
good  idea.  But,  it  cannot  be  done  without  an  understanding  of  the  uncertainty  in  the  density  functional.  We 
suggest  that  one  considers  using  the  semi-empirical  density  functionals  that  the  chemists  have  constructed 
by  fitting  a  guessed  at  form  with  parameters  to  the  results  of  configuration  interaction  computations  for  small 
systems  [4,  7].  This  can  be  taken  as  a  standard  for  comparisons.  Then  the  various  physically  motivated  func¬ 
tionals  would  deviate  from  that  standard,  and  one  could  study  the  propagation  of  those  deviations  through 
the  Car-Parrinello  computation.  It  has  become  apparent  to  solid-state  physicists  that  no  functional  or  com¬ 
putational  method  that  does  not  treat  exchange  exactly  is  to  be  taken  seriously.  This  is  an  excellent  problem 
for  uncertainty  analysis  in  molecular  modelling. 

4.  Molecular  Systems. 

Interactions  in  molecular  systems  can  be  divided  into  intra-molecular  and  inter-molecular  ones. 
The  latter  are  often  taken  to  be  atom-pair  interactions.  The  intra-molecular  interactions,  i.e.,  the  interactions 
between  the  atoms  of  one  molecule,  are  determined  by  chemical  bonds  and  are  therefore  not  only  strong 
compared  with  inter-molecular  interactions  (between  atoms  of  different  molecules)  but  also  include  oriental 
dependence.  A  brief  description  of  intra-molecular  degrees  of  freedom  and  interactions  is  as  follows. 

First  of  all,  the  chemical  bonds  can  stretch.  The  interaction  associated  with  this  degree  of  freedom 
is  usually  described  in  the  form  of  a  harmonic  potential  for  the  bond  length  / 

Vs„e,ch«)  =  (22) 


where  l0  is  the  equilibrium  bond  length. 

Now  consider  three  atoms  bonded  in  a  chain-like  configuration  A-B-C.  This  chain  is  characterised 
by  a  bending  or  valence  angle  (p  which  varies  around  an  equilibrium  value  cpQ  and  the  potential  is  described 
in  terms  of  a  cosine 


V valence^  =  -CCp[C0S((p -  q>0)  +  COS((p  +  (pQ)] 


(23) 
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where  the  equivalence  of  the  angles  cp0  and  ~<Po  taken  into  account. 

Finally  there  is  an  interaction  associated  with  chain  configurations  of  four  atoms  A-B-C-D.  The 
plane  through  the  first  three  atoms,  A,  B ,  and  C  does  not  in  general  coincide  with  that  through  B,  C,  and  D. 
The  torsion  interaction  is  similar  to  the  bend  interaction,  but  the  angle,  denoted  by  0,  is  now  between  these 
two  planes 


v  torsion^)  =  -ar[cos(e-e0)  +  cos(e  +  60)]  (24) 

Characteristic  vibrations  associated  with  the  different  degrees  of  freedom  distinguished  here  can  be  derived 
from  the  harmonic  interactions  -  the  frequencies  vary  as  the  square  root  of  the  a-coefficients.  In  general  the 
bond  length  vibrations  are  the  most  rapid,  followed  by  the  bending  vibrations.  In  order  for  molecular 
dynamics  simulations  to  be  accurate,  the  time  step  for  integration  should  be  chosen  smaller  than  the  fastest 
degree  of  freedom.  But  this  degree  of  freedom  will  vibrate  with  a  small  amplitude,  because  of  the  strong 
potential,  most  of  the  computational  time  in  molecular  dynamics  simulations  will  be  consumed  by  those 
parts  of  the  motion  which  are  not  expected  to  contribute  strongly  to  the  physical  properties  of  the  system. 
Moreover,  if  there  is  a  clear  separation  between  the  time-scales  of  the  various  degrees  of  freedom  of  the  sys¬ 
tem,  energy  transfer  between  the  fast  and  slow  modes  is  extremely  slow,  so  it  is  difficult,  if  not  impossible, 
to  reach  equilibrium  within  a  reasonable  amount  of  time.  It  is  in  such  instances  of  molecular  dynamics  that 
graph-theoretic  methods  for  the  identification  of  almost  invariant  sets  can  simplify  and  accelerate  conver¬ 
gence  to  equilibrium  [8]. 

4.1  Direct  Approach  for  the  Nitrogen  Molecule. 

A  good  base  case  for  investigating  methods  of  uncertainty  analysis  is  molecular  dynamics  of  the 
nitrogen  molecule  N2  in  the  context  of  the  rigid  molecule  approximation.  In  this  approximation  molecules 
are  treated  as  rigid  bodies  whose  motion  consists  of  translations  of  the  centre  of  mass  and  rotations  about 
this  point.  The  forces  acting  between  two  rigid  molecules  are  usually  composed  of  atomic  pair  interactions 
between  atoms  belonging  to  different  molecules.  The  total  force  acting  on  a  molecule  determines  the  transla¬ 
tional  motion  and  the  torque  determines  the  rotational  motion.  The  N2  molecule  consists  of  two  nitrogen 
atoms,  each  of  mass  14  atomic  mass  units  (amu)  and  whose  separation  d  is  kept  fixed  in  the  rigid  approxima¬ 
tion.  The  coordinates  of  the  molecule  are  the  three  coordinates  of  the  centre  of  mass  and  the  two  coordinates 
defining  its  orientation.  The  latter  can  be  polar  angles  but  here  we  shall  characterise  the  orientation  of  the 
molecule  by  a  unit  normal  pointing  ft  from  one  atom  to  the  other. 

The  motion  of  the  centre  of  mass  of  the  molecule  is  determined  by  the  total  force  Ftot  acting  on  the 
particular  molecule.  This  force  is  the  sum  of  all  the  forces  between  each  of  the  two  atoms  in  the  molecule 
and  in  the  remaining  molecules.  The  atomic  forces  can  be  modelled  by  a  Lennard-Jones  interaction  with  the 
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appropriate  atomic  nitrogen  parameters  a  =  3.72  Angstroms  and  z/kB  =  37.3  Kelvin  [9]  where  kB  is  the 
Boltzmann  constant.  The  equation  of  motion  for  the  centre  of  mass  Rcm  is  then 

*cm  =  Ftot  (25) 

which  can  be  solved  in  exactly  the  same  way  as  in  an  ordinary  molecular  dynamics  simulation.  The  motion 
of  the  orientation  vector  ft  is  determined  by  the  torque  x  with  respect  to  the  centre  of  the  molecule  which  is 
given  in  terms  of  the  forces  and  F ^  acting  on  atoms  1  and  2  respectively. 

T  =  (d/2)hx(Fm-Fm).  (26) 

The  torque  changes  the  angular  momentum  L  of  the  molecule.  This  is  equal  to  /co,  where  I  is  the  moment  of 

inertia  equal  to  md2  and  co  is  the  angular  frequency  vector  whose  norm  is  the  angular  frequency  and  whose 
direction  is  the  axis  around  which  the  rotation  takes  place.  Note  that  x  is  not  necessarily  parallel  to  CD.  The 
equation  of  motion  for  the  angular  momentum  is 

x  =  /(b  (27) 

and  we  note  that  the  angular  frequency  CD  is  in  turn  related  to  the  time  derivative  of  the  direction  vector  ft  as 

in  =  (OXh  (28) 

at 

Combining  equation  (27)  and  (28),  the  molecular  dynamics  equation  is 

d 1  2 

-~h  =  co  x  (co  x  «)  +  x  x  h/I  =  -  CD  ft  +  x x  h/1  (29) 

dr 

and  this  equation  of  motion  leaves  the  norm  of  the  direction  vector  h  as  it  should. 

For  general  molecules,  there  will  exist  an  extra  degree  of  freedom:  the  angle  of  rotation  around  a 
molecular  axis  -  the  third  Euler  angle,  which  is  denoted  as  y.  The  straightforward  procedure  for  solving  the 
equations  of  motion  is  to  calculate  the  principal  angular  velocity  co  in  terms  of  the  time  derivatives  of  the 
Euler  angles.  The  Euler  equation  of  motion  gives  the  rate  of  change  in  co  in  terms  of  the  torque.  The  time 
derivatives  of  the  Euler  angles  can  be  found  again  from  co,  and  these  can  be  used  to  calculate  the  new  atomic 
positions.  There  is  however,  a  problem  when  the  Euler  angle  0  is  zero,  as  in  that  case  the  transformation 
from  co  to  the  time  derivatives  of  the  Euler  angles  becomes  singular.  The  most  efficient  method  is  to  use  the 
quaternion  representation  which  avoids  the  instability  resulting  from  this  singularity.  In  the  quaternion 
method,  the  orientation  of  the  molecule  is  represented  in  terms  of  a  four-dimensional  unit  vector  rather  than 
three  Euler  angles  [11]. 
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4.2  The  CS2  Trimer. 

A  more  complicated  example  is  the  trimer  molecule  CS2  [12].  The  linear  geometry  of  this  molecule 
is  in  principle  imposed  automatically  by  the  correct  bond  constraints  between  the  three  pairs  of  atoms.  How¬ 
ever,  the  motion  of  this  molecule  is  described  by  two  positional  degrees  of  freedom:  two  to  define  the  orien¬ 
tation  of  the  molecule  and  three  for  the  center  of  mass  position.  The  three  atoms  without  constraints  have 
three  degrees  of  freedom  and  if  three  of  these  are  eliminated  using  the  bond  constraints,  we  are  left  with  six 
degrees  of  freedom  instead  of  the  five  required.  A  procedure  is  to  fix  only  the  distance  between  the  two  sul¬ 
phur  atoms 


\rtfi)~rs(2)\2  =  ^  (30) 

and  to  fix  the  position  of  the  C-atom  by  a  linear  vector  constraint 

(rsn)  +  r<4\))/2~rc  -  0  (31) 

adding  up  to  the  four  constraints  required.  For  a  molecule,  in  general  a  number  of  atoms  forming  a  backbone 
set  is  identified  and  these  are  fixed  by  bond  constraints  (the  two  sulphur  atoms  in  this  example).  In  the  case 
of  a  planar  structure  we  take  three  non-colinear  atoms  as  a  backbone.  These  atoms  satisfy  three  bond  con¬ 
straints  and  the  remaining  atoms  are  constrained  linearly.  In  a  three-dimensional  molecular  structure,  four 
backbone  atoms  are  subject  to  six  bond  constraints  and  the  remaining  ones  to  a  linear  vector  constraint  each. 
In  the  constraint  procedure,  the  degrees  of  freedom  of  the  non-backbone  atoms  are  eliminated  so  that  the 
forces  they  feel  are  transferred  to  the  backbone.  This  elimination  is  always  possible  for  linear  constraints 
such  as  those  obeyed  by  the  non-backbone  atoms. 

The  equations  of  molecular  dynamics  for  this  system  may  be  written  down  now  by  introducing 
Lagrange  multipliers  X  and  [i  for  the  bond  and  linear  vector  constraints  respectively. 


2 

drSl"  17  VW  .  |X 

ms— 2  =  ^1  -  -  rs(2))  - 1 


(32) 


2 

dr 


Sa) 

ms — 2  ”  ^2  +  ~~  rc<2)) 

dt 


sfl)~  sf2”  -  2 


(33) 


m 


2 

drC 


(34) 


The  linear  constraint  of  equation  (3 1 )  is  now  differentiated  twice  with  respect  to  time  and  using  the  equa¬ 
tions  of  motion 
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Fc^  =  itSF  1+V^ 


(35) 


is  obtained.  We  can  thus  eliminate  the  multiplier  |i  and  obtain  the  equations  of  motion  for  the  S-atoms 


df 


(  mr\  mc  mc 

('-mr, 


(36) 


d  r 


df 


(  mr\  mr  mr 


(37) 


where  M-2ms  +  mc.  It  is  not  possible  to  eliminate  the  multiplier  X.  Therefore,  the  molecular  dynamics 

equations  for  CS2  involve  integrating  equations  (36)  and  (37)  subject  to  the  constraint  of  equation  (30).  This 
is  achieved  by  the  modified  Verlet  algorithm  of  Section  3.1  with  the  atomic  forces  F}  and  F2  being  appropri¬ 
ately  negative  gradients  of  the  interatomic  potential  [13]. 

Bibliography. 


1.  Ashcroft,  N.  W.  and  N.  D.  Mermin,  1976,  Introduction  to  Solid  State  Physics ,  Holt,  Reinhart,  and 
Wisnton,  New  York. 

2.  Allen,  M.  P.  and  D.  J.  Tildesley,  1989,  Computer  Simulation  of  Liquids,  Oxford  University  Press. 

3.  Rountree,  C.  L.,  Kalia,  R.  K.,  Lidorikis,  E.,  Nakano,  A.,  van  Brutzel,  L.,  and  P.  Vashishta,  2002, 
Annual  Reviews  of  Materials  Research ,  377-400. 

4.  Leach,  A.  R.,  2001,  Molecular  Modelling  -  Principles  and  Applications ,  Second  Edition,  Prentice 
Hall. 

5.  Kohn,  W.,  1999,  Reviews  of  Modern  Physics ,  71,  1253-1266. 

6.  Car,  R.  and  M.  Parrinello,  1985,  Phys .  Rev.  Lett.,  55,  2471. 

7.  Pople,  J.  A.,  1999,  Reviews  of  Modern  Physics,  1 1 ,  1267-1274. 

8.  Dellnitz,  M.  and  R.  Preis,  2002,  Congestion  and  Almost  Invariant  Sets  in  Dynamical  Systems,  and  ref¬ 
erences  therein. 

9.  Cheung,  P.  S.  Y.  and  J.  G.  Powles,  1975,  Molecular  Physics,  30,  921-949. 

10.  Evans,  D.  J.,  1977,  Molecular  Physics,  34,  317-325. 

1 1 .  Evans,  D.  J.  and  S.  Murad,  1977,  Molecular  Physics,  34,  327-331 . 

12.  Ciccotti,  G.,  Ferrario,  M.,  and  J.  P.  Ryckaert,  1982,  Molecular  Physics,  47,  1253-1264. 

13.  Goldstein,  H.,  1980,  Classical  Mechanics,  Addison-Wesley,  Reading. 


Symmetry-Breaking  and  Uncertainty  Propagation  in  a  Reduced  Order  Thermo-acoustic  Model 


Gregory  Hagen  Andrzej  Banaszuk 


Abstract — We  present  a  thermo-acoustic  model  on  a  cylin¬ 
drical,  or  annular,  geometry,  capable  of  modeling  instabilities 
of  tangential  acoustic  modes.  The  model  accounts  for  non- 
uniform  density,  damping,  rotational  flow,  and  heat-release 
coupling.  It  is  shown  that  deliberately  introducing  spatial 
variations  in  some  quantities  has  a  similar  effect  to  adding 
damping  to  the  system.  The  effects  of  these  symmetry-breaking 
conceptes  are  evaluated  on  the  model  through  linear  analysis 
and  the  net  amount  of  additional  damping  is  computed.  We 
show  how  various  symmetry-breaking  concepts  are  robust  with 
respect  to  the  uncertainty  in  the  model  parameters  and  we 
examine  propagation  of  uncertainty  with  respect  to  a  recently 
defined  measure  of  uncertainty. 


A 

a 

e 

K 

P 

Q 

r,e 

R 

t 

u 

Ur  j  Uq 

v2{w) 

X 

Y 

7 

P 

<P 

$ 

£ 

(■)>(•)' 

(:) 

(’)Oj  (')m 


Nomenclature 

=  state  space  matrix 

—  speed  of  sound 

=  internal  energy  per  unit  volume 
=  heat-release  gains 
=  pressure 

=  volumetric  heat-release 
=  radial,  tangential  (space)  variables 
=  mean  radius 
=  time 

=  velocity  vector 
=  radial,  tangential  velocities 

—  uncertainty  metric 

=  state  vector  in  Galerkin  truncated  model 
=  acoustic  boundary  admittance 

—  ratio  of  specific  heats  (=  1.4) 

=  density 

—  tangential  (Fourier)  basis  function 
=  radial  basis  function 

=  damping  constant 
=  mean,  perturbation  quantities 
=  normalized  quantities 
=  spatially  mean  and  m-periodic  quantities 


I.  Introduction 

Thermo-acoustic  instabilities  in  gas  turbine  and  rocket 
engines  develop  when  acoustic  waves  in  combustors  couple 
with  an  unsteady  heat-release  field  in  a  positive  feedback 
loop.  For  a  summary  of  active  control  of  thermo-acoustic 
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instabilities  see  [3].  Thermo-acoustic  modeling  and  control 
is  well-studied  for  axially  extended  combustion  chambers, 
as  in  [5],  [7],  [14],  [10],  where  the  acoustic  to  heat- 
release  coupling  is  dominated  by  longitudinal  acoustic 
modes.  However,  comparably  less  attention  has  focussed 
on  thermo-acoustic  modeling  in  combustion  chambers  with 
annular,  or  cylindrical  geometries.  In  the  present  work,  we 
develop  a  low-order  thermo-acoustic  model  on  a  circular 
geometry,  somewhat  similar  to  that  of  [1]. 

Often  thermo-acoustic  instabilities  are  dominated  by  a 
few  natural  acoustic  modes,  which  can  accurately  be  mod¬ 
eled  with  a  low-dimensional  model.  Accurate  heat-release 
models  are  difficult  and  time-consuming  to  implement.  A 
low-order  thermo- acoustic  model,  properly  calibrated  with 
acoustic  data,  can  provide  insight  into  the  possibly  dele¬ 
terious  acoustic-heat-release  coupling  and  may  provide  a 
platform  for  fast  evaluation  of  preliminary  design  concepts. 

In  this  work  we  define  symmetry-breaking  as  the  de¬ 
liberate  introduction  of  spatial  variations  in  the  system  pa¬ 
rameters  in  order  to  change  the  stability  properties.  Recent 
work  has  focussed  on  analysis  of  heterogeneous  distributed 
systems  [6], [9], [11].  Symmetry-breaking  is  commonly  ref- 
ered  to  as  mistuning  in  the  literature  regarding  the  dynamics 
of  arrays  of  turbine  blades  on  a  disk.  Studies  of  stability 
properties  of  turbine  blade  flutter  through  the  introduction 
of  spatial  nonuniformities  has  appeared  in  [2],  [16].  Optimal 
mistuning  in  arrays  of  bladed  disks  has  appeared  in  [15], 
[17].  A  study  of  the  effects  of  asymmetry  on  compressor 
stall  inception  has  appeared  in  [8]. 

As  in  the  case  of  mistuning  in  arrays  of  bladed  disks  in 
turbines,  this  form  of  passive  control  is  often  more  feasible 
than  implementing  an  active  control  scheme.  This  may 
also  be  true  for  the  case  in  combustion  chambers,  where 
high  temperatures  prohibit  adequate  sensing  and  may  dam¬ 
age  the  actuators  required  for  active  control.  Furthermore, 
symmetry-breaking  can  be  a  more  cost-effective  means  of 
stability  enhancement. 

In  the  present  work  we  develop  a  thermo-acoustic  model 
capable  of  describing  oscillations  and  passive  control  of 
tangential  modes  using  symmetry-breaking  concepts.  The 
model  is  derived  in  section  II,  where  we  describe  the  param¬ 
eters  available  for  symmetry-breaking  for  stability  enhance¬ 
ment.  The  low-order  model  is  constructed  by  Galerkin  pro¬ 
jection  and  truncation.  Linear  analysis  of  stabilization  via 
symmetry-breaking  is  discussed  in  section  III.  Parametric 
uncertainties  are  modeled  as  spatially  varying  components 
within  the  model.  The  propagation  of  these  uncertainties 


to  the  amplitude  of  pressure  mode  oscillations  is  described 
in  terms  of  a  recently  defined  measure  of  uncertainty  [12], 
which  is  discussed  in  section  IV. 


II.  Model  Development 


We  consider  a  combustor  that  has  an  annular  shape  and 
includes  a  circumferential  array  of  bluff  body  flame  holders, 
similar  to  that  presented  in  [1].  Flameholders  extend  radially 
from  inner  to  outer  diameter  of  the  annular  combustor.  A 
cut  along  a  contant  radius  surface  is  shown  in  Figure  1. 
The  thermo-acoustic  system  we  consider  is  the  linearized 


Fig.  1.  A  radial  cross-section  of  a  bluff  body  flameholder  array  in  the 
tangential  direction. 


momentum  and  pressure  dynamics  as  derived  by  Culick  [4]. 
We  now  list  the  assumptions  used  throughout  the  rest  of  this 
work: 

•  Viscous  effects  are  negligible. 

•  All  quantities  are  periodic  in  0. 

•  The  acoustic  medium  is  air  and  it  acts  as  an  ideal  gas, 
with  constant  ratio  of  specific  heats,  7. 

•  The  flow  is  non-conducting  (i.e.  the  other  transport 
terms  appearing  in  [4]  are  negligible). 

•  The  steady  density  p(9),  tangential  velocity  u$(0), 
and  boundary  admittance  Y(9)  are  possibly  spatially 
varying. 

•  All  other  steady,  mean  quantities  are  spatially  uniform. 

•  The  radial  velocity  perturbation  ulr  is  negligible  on  the 
interior  of  the  domain. 

•  Variations  in  the  axial  direction  are  negligible. 

•  The  heat-release  and  damping  feedback  are  time- 
invariant. 

•  The  heat-release  feedback  is  a  scalar  operation  that 
depends  only  on  the  tangential  acoustic  velocity. 


We  start  with  the  inviscid  transport  equations  written  in 
cylindrical  coordinates, 

dp 


dt+v-(»p)  = 

<9u 

P-37+Pu'Vu  = 


0 


dt 
de 

PTt+pn-Ve 


—  —S/p 
=  -pS/  •  u  +  q. 


(1) 

(2) 

(3) 


By  substitution  of  the  ideal  gas  relation  e  =  ^  an^  (1) 

we  obtain 


(7-1) 


de 

p-+pu.Ve 
dp  p  dp  —  p  _ 

at--pa!  +  a7p--fu'Vp 

+  -V  •  {up)  +  u  •  Vp  -  -u  •  Vp 

Ot  p  p 


dp 

dt 


-f-  pSJ  ■  u  +  u  •  S/p. 


(4) 


Multiplying  (3)  by  (7  • 
dp 


1)  and  applying  (4)  results  in 

1  )<F  (5) 


•tz  +  u  •  S/p  +  7 pSZ  .  u  =  (7 
ot 

The  system  (2,5)  is  the  thermo-acoustic  model  that  we 
focus  on.  Equations  (2,5)  are  linearized  about  the  steady 
pressure  p,  density  p(9),  and  tangential  velocity  u$(9) 
where  p  is  spatially  uniform,  and  p(0)  and  ue(0)  are 
possibly  spatially  varying  in  the  tangential  direction.  For 
notational  convenience,  we  drop  the  arguments  (0)  in  the 
following  equations.  The  linearized  dynamics  are 


dv'r  _  _puedK  2pue 
dt  r  dO  r 


u„  - 


1  <y 

dr 


_du'e  __  puo  du'e  pv!e  due  P^e  ,  1  dp' 
^  dt  r  dO  r  dO  r  Ur  r  dO 


(6) 


(7) 


dp'  *ip  du'e  7 p'  due  7 p  dru'r 
dt  ^  r  dO  ^  r  dO  r  dr 

__3g +  ,8) 

For  simplicity,  we  assume  that  u'r  is  negligible  on  the 
interior  of  the  domain,  however,  we  retain  that  u'r  couples 
with  the  pressure  on  the  outer  boundary  through  a  static 
boundary  admittance  relations  (see  e.g.  [13]) 

<  =  Y(0)p,  (9) 

where  the  argument  (0)  is  included  to  indicate  possible 
spatial  variations  in  the  boundary  admittance.  With  the 
exception  of  the  admittance  relation  (9),  we  assume  that 
the  contribution  of  the  radial  acoustic  velocity  coupling  to 
the  dynamics  of  the  thermo-acoustic  system  is  negligible, 
hence  we  will  ignore  (6).  Note  that  further  generalizations 
of  this  model  can  retain  these  dynamics.  Also,  note  that  a 
more  general  linear  dynamic  admittance  condition  could  be 
easily  handled  by  frequency  domain  techniques. 


( 


We  will  assume  some  generic,  normalized  mode  shape 
in  the  radial  direction,  tp(r)  which  may  be  some  weighted 
combination  of  Bessel  functions,  depending  on  the  radial 
boundary  conditions  [4]  ,  We  assume  that  each  quantity  is 
written 


We  project  the  sysetm  on  to  this  radial  mode.  Integrating 
the  term  with  vlr  in  (8)  results  in 


ty{R<2)R2v!T{R<2) 

iKWalWa) 

4>(R2)2R2Yp'(0,t),  (10) 


where  we  have  applied  (9)  and  the  assumption  that  u'r  =  0 
on  the  interior.  Equations  (7)  and  (8)  can  be  simplified  to 


_dufQ  puo  du'e  pule  due  1  dp' 
P~8t  =  ~~W~86  ~  -R~de  -  R~de 


dt 


R  86  W 


1)9', 


(12) 


where  R  is  the  mean  radius.  For  notational  convenienc  we 
define  the  normalized  quantities 


u  -  Ug, 


(13) 


£  =  ypip(R2)2R2Y. 
The  speed  of  sound  is  given  by 


(14) 


Spatially  Varying  Parameters 

We  expand  all  of  the  quantities  in  terms  of  Fourier  modes 
{0}n  where 

So  the  pressure  is  expanded  as 

P=y2Pn'  Pn-(pAn)  (19) 

n 

where  (-,*)  is  the  standard  inner  product  on  L2(0,  27t). 
The  coefficients  (*)n  are  defined  similarly  with  respect  to 
all  of  the  other  parameters  in  the  model.  We  assume  that 
the  spatially  varying  parameters  have  period  m  around  the 
circle  in  addtion  to  some  mean  quantity.  Thus,  we  define 


Ug(0) 

Uoo  +  ^6m  [0m (0)  +  </)-m(9j\ 

(20) 

m 

'=  $0  +  [0m(^)  +  4>-m{9j] 

(21) 

a2(0) 

:=  a0  +  am  [0m (^)  +  0-m(0)]  • 

(22) 

We  assume  that  the  heat-release  feedback  is  given  by 
q(u)  =  Ku.  The  case  where  q(u)  is  nonlinear  will  be 
examined  in  section  IV.  The  spatially  varying  heat-release 
gain  is  given  by 

K(0)  :=  K0  +  Km  [<f>m(9)  +  0_m(0)] .  (23) 


Galerkin  Projection 

We  are  interested  in  the  stability  of  the  first  tangential 
mode,  and  therefore  we  apply  a  Galerkin  projection  of  the 
system  (17,  18)  on  to  the  first  Fourier  modes  +1,-1.  We 
have  [0m  +  0_m]  [01  +  0—1  ]  =  0771+1  +  0m— 1  +  0-m+i  + 
0_m-i.  Since  we  are  focussing  on  the  stability  of  the  first 
mode  only,  we  take  m  —  1  in  the  (20-23).  This  produces 
additional  cross-coupling  between  mode  +1  and  mode  —  1. 


a 


2 


yp 

~P  * 


(15) 


Because  ~p  is  possibly  spatially  varying,  the  speed  of  sound 
is  also  spatially  varying  and  we  write  a2  =  a2  (6).  Note  that 
p  is  now  dimensionless  and  q  has  dimensions  of  [^—].  We 
scale  9  as 

0  ->  Re,  (16) 


The  Galerking  truncation  of  the  system  (17,  18)  is 


x  =  Ax. 


(24) 


where  x  =  [u\,  ti_i,pi,p_i]T  and 


A  = 


-tug  o  -lUg  2  ~ia%  2 
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so  9  now  has  dimensions  of  [length].  We  apply  (13-16)  to 
(11,12)  to  result  in  the  system 
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(17) 
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for  9  e  [0,  27t).  Equation  (17)  has  units  of  and 

equation  (18)  has  units  of  [g— ]. 


(25) 

HI.  Linear  Analysis  of  Symmetry-Breaking 

In  this  section  we  present  some  computational  results  that 
evaluate  effect  of  various  symmetry-braking  concepts  on  the 
eigenvalues  of  (25).  In  the  following  computational  results, 
nominal  values,  £0  =  0.2,  ue0  —  0.05,  ag  —  1,Kq  =  0.5 
were  used.  In  the  following  analysis,  variations  in  the  mean 
tangential  flow  had  little  effect  on  the  stability  of  the  system. 


Damping  Variation 

Figure  2  shows  the  maximum  real  parts  of  the  stable 
eigenvalues  of  A  resulting  from  differenct  combinations  of 
£o  and  £2.  A  value  of  zero  in  the  figure  indicates  that  all  of 
the  eigenvalues  are  in  the  right  half  of  the  complex  plane. 
For  a  fixed  value  of  uniform  damping  £0,  increasing  £2 
decreases  the  overall  damping  of  the  system.  Therefore, 
uniform  damping  is  more  effective  than  spatially  varying 
damping. 


maximum  real  part  of  eigenvalues 


Fig.  2.  Maximum  real  parts  of  stable  eigenvalues  for  varying  fo  and  £2- 
A  value  of  zero  indicates  that  the  linear  system  is  unstable. 


Speed  of  Sound  Variation 

Variation  in  the  speed  of  sound  is  equivalent  to  changing 
the  local  resonance  frequency.  The  results  summarized  in 
[2]  show  that  such  variations  can  be  optimized  to  increase 
damping.  Figure  3  shows  the  maximum  real  parts  of  the  sta¬ 
ble  eigenvalues  of  A  resulting  from  differenct  combinations 
of  al  and  a|.  The  figure  shows  that  increasing  the  spatial 
variation  of  the  speed  of  sound  has  a  comparable  effect  to 
increasing  £0  (Compare  with  Figure  2).  We  see  that  changes 
in  al  on  the  order  of  10%  of  have  a  significant  effect 
on  the  stability  of  A. 

Heat-Release  Gain  Variation 

Spatial  variations  in  the  heat-release  coupling  is  an  exam¬ 
ple  of  dynamic  symmetry-breaking  in  the  system.  In  other 
words,  it  may  be  possible  to  change  the  dynamic  heat- 
release  coupling,  while  maintaining  mean  quantities,  such 
as  temperature,  density,  etc.  Figure  4  shows  a  root  locus 
plot  of  the  eigenvalues  resulting  from  variations  in  K2  and 
Kq.  We  see  that  increasing  the  spatial  variation  of  the  heat- 
release  gain  moves  two  eigenvalues  to  the  left  and  moves 
the  other  two  eigenvalues  to  the  right.  A  contour  plot  similar 
to  Figure  3  can  be  obtained  for  general  combinations  of  K2 
for  a  constant  Kq . 


Fig.  3.  Maximum  real  parts  of  stable  eigenvalues  for  varying  and  a%. 
A  value  of  zero  indicates  that  the  linear  system  is  unstable. 
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Fig.  4.  Root  locus  plot  showing  eigenvalues  for  varying  values  of 
K2.  Eigenvalues  corresponding  with  K2  =  0  are  shown  in  blue,  and 
eigenvalues  corresponding  with  K2  =  0.5  are  shown  in  blue. 


IV.  Uncertainty  Propagation 

We  analyze  uncertainty  propagation  with  respect  to  the 
uncertainty  measure,  denoted  by  u(-),  recently  defined  in 
[12].  For  the  present  study,  the  output  observable  of  interest 
is  the  amplitude  of  oscillation  of  the  first  tangential  pres¬ 
sure  mode.  We  denote  by  vj(p)  the  resulting  probability 
distribution  of  the  pressure  amplitude  which  results  from 
a  given  distribution  of  uncertain  model  parameters.  The 
uncertainty  measure  is  therefore  denoted  by  v(m).  For  the 
case  of  no  uncertainty  (all  model  parameters  are  exactly 
known)  the  resulting  probability  distribution  will  be  a  dirac- 
delta  function  and  the  resulting  uncertainty  metric  will  be 
v(vj)  —  0. 

As  an  example,  we  examine  uncertainty  propagation  with 
respect  to  the  measure  vj(p),  that  results  from  uncertain 
spatially  varying  sound  speed,  a\  with  different,  fixed  and 
certain,  values  of  mean  heat-release  gain  Kq.  Similar  results 
can  be  obtained  for  all  of  the  other  combinations  of  model 
parameters.  Figure  5  shows  the  maximum  real  parts  of  the 
stable  eigenvalues  of  A,  with  nominal  values  a§  =  l,£o  = 


Histograms  (bars)  and  Cumulative  Distribution  Functions  (lines) 


0.2.  We  examine  cases  with  Ko  varying  from  0.1  to  0.35 
where  a 2  is  uncertain  and  may  be  up  to  20%  of  a2. 


We  compute  the  resulting  uncertainty  metric  as  defined 
in  [12], 


V2(rv)  =  min \w  -  Sz\ 2, 


where  |  •  |2  denotes  the  L2- norm  on  the  space  M  of 
probabilistic  measure  on  M.  It  is  convenient  to  use  the  cumu¬ 
lative  distribution  function  F to  compute  the  uncertainty 
measure.  In  this  case,  since  FVJ(p)  is  continuous,  it  can  be 
shown  that  the  uncertainty  measure  can  be  computed  by 


We  present  the  results  of  Monte-Carlo  simulations  of 
the  PDE  system  (17,  18)  where  the  heat-release  coupling 
includes  a  saturation  effect.  From  the  analysis  of  the  trun¬ 
cated  model  and  Figure  5,  we  see  that  as  #0  increases,  the 
system  approaches  the  linear  stability  boundary,  and  hence, 
the  amplitude  should  increase. 


maximum  real  part  of  eigenvalues 


Ko 


Fig.  5.  Maximum  real  parts  of  stable  eigenvalues.  A  value  of  zero 
indicates  that  the  linear  system  is  unstable. 


The  resulting  distributions  w{p)  (shown  as  bars)  and 
the  respective  normalized  cumulative  distribution  functions 
Fzu(p)  (shown  as  lines)  of  the  pressure  amplitude  are  shown 
in  Figure  6,  for  the  nominal  heat-release  gain  Ko  ranging 
from  0.1  to  0.35.  For  small  vales  of  K0  the  system  is 
linearly  stable  and  therefore,  all  of  the  resulting  amplitudes 
are  zero.  Therefore,  the  distribution  of  amplitudes  for  these 
cases  is  a  delta  function  at  zero.  Because  of  this,  the 
uncertainty  propagation  in  terms  of  the  metric  (26)  is  also 
zero.  As  Ko  increases,  the  system  approaches  the  linear 
stability  limit,  and  the  resulting  distribution  of  amplitudes 
has  a  wider  range.  Therefore  the  uncertainty  propagation 
corresponding  with  these  cases  is  larger.  When  Ko  increases 
beyond  the  stability  limit,  the  system  enters  a  limit-cycle 


Fig.  6.  Histograms  and  cumulative  distribution  functions  of  resulting 
amplitudes  of  oscillation  for  varying  Kq  and  uncertainty  in  a|. 


where  the  amplitude  of  oscillation  is  determined  by  the 
saturation  nonlinearity  and  the  amount  of  damping. 

Figure  7  shows  the  values  of  the  uncertainty  metric  based 
on  the  observable  of  pressure  amplitude  corresponding  with 
different  values  of  Ko  and  uncertainty  in  a 2  up  to  20% 
of  a2.  We  see  that  the  lowest  amount  of  uncertainty  in 
the  amplitude  correspond  with  the  highly  stable  cases.  The 
uncertainty  measure  increases  as  the  system  approaches 
the  stability  boundary.  The  highest  amount  of  uncertainty 
occurs  near  the  linear  stability  boundary  of  A.  However, 
for  high  values  of  Ko ,  the  uncertainty  propagation  of  a2  to 
the  pressure  amplitude  is  actually  decreased.  We  conjecture 
that  the  increased  stability  of  the  limit  cycle  explains  the 
reduced  uncertainty  measure  for  the  limit-cycling  cases. 
An  analytic  study  to  investigate  the  relationship  between  a 
measure  of  hyperbolicity  of  the  dynamics  and  the  measure 
of  uncetrainty  is  currently  in  progress. 


Fig.  7.  Uncertainty  metric  v{zv)  versus  Kq. 


V.  Conclusion 

We  have  presented  a  simple  one-dimensional  thermo¬ 
acoustic  model  that  is  adequate  for  modeling  oscillations 
of  tangential  acoustic  modes.  The  model  has  been  derived 
allowing  for  some  spatially- varying  parameters.  We  have 
investigated  the  effects  of  symmetry-breaking,  or  the  delib¬ 
erate  introduction  of  spatial-variations  in  system  parameters, 
on  the  system’s  linear  stability  properties.  It  was  found  that 
spatial  variations  in  the  steady  tangential  flow  and  damping 
do  not  enhance  stability.  Variations  in  the  local  speed 
of  sound  and  the  heat-release  coupling  gains  significantly 
affects  system  damping.  An  example  of  uncertainty  propa¬ 
gation  through  this  system  was  presented  in  the  context  of 
spatially  varying,  uncertain,  speed  of  sound,  with  varying 
levels  of  mean  heat-release  coupling.  In  this  example,  the 
uncertainty  propagation  given  in  terms  of  an  uncertainty 
measure  depending  on  the  resulting  distribution  of  the  pres¬ 
sure  amplitude,  was  the  lowest  for  the  limit-cycling  case. 
An  analytic  study  to  investigate  the  relationship  between  a 
measure  of  hyperbolicity  of  the  dynamics  and  the  measure 
of  uncetrainty  is  currently  in  progress. 
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8:00  am:  Introduction  to  UTRC,  Dr.  Jean  Colpin,  Director,  UTRC 

8:30  am:  Opening  Remarks,  Dr.  Carey  Schwartz,  Program  Manager,  DARPA 

9:00  am:  Analytical  Systems  Engineering,  Dr.  Clas  Jacobson,  Director,  Systems  Dept,  UTRC 

9:30  am:  Design  of  Robust  Dynamical  Systems:  Key  Principles  and  Workshop  Overview 

Andrzej  Banaszuk,  Pratt  &  Whitney  Program  Office,  UTRC 

10:00  am  Break 

Uncertainty  Propagation  in  Dynamical  Systems 

10:30  am:  Uncertainty  in  Analysis  &  Design:  Dynamical  Systems  Perspective 

Igor  Mezic,  Professor,  Dept.  Mechanical  &  Environmental  Engineering,  UCSB 

11:30  am:  Uncertainty  and  networks  of  dynamical  systems 

Jerrold  Marsden,  Professor,  Control  &  Dynamical  Systems,  Caltech 

12:30  pm:  Working  Lunch 

1:30  pm:  Uncertainty  Propagation  in  networks  of  Dynamical  Systems 

Tamas  Kalmar-Nagy,  Mihai  Huzmezan,  UTRC 

Decompositions  of  Large  Networks  of  Dynamical  Systems 

2:00  pm:  Analysis  and  Computational  Methods  for  Dynamical  Systems  on  Graphs 

Sanjay  Lall,  Stanford  and  Matthew  West,  UC  Davis 

3:00  pm:  On  The  Approximation  of  Complicated  Dynamics,  Michael  Dellnitz,  Department 
Chair,  Applied  Mathematics,  University  of  Paderbom 


4:00  pm:  Break 


4:30  pm:  Graph  Decomposition  Methods 

Robert  LaBarre,  Subbarao  Varigonda,  UTRC 

5:00  pm  Propagation  of  Invariant  Measures  in  Interconnected  Dynamical  Systems 

Leonid  Bunimovich,  Professor,  Georgia  Institute  of  Technology 


6:00  pm:  Dinner 

Friday,  January  16 

7:30  am:  Continental  Breakfast 

Iterative  Methods  for  Dynamical  Systems 

8:00  am:  Spectral  Balance:  Frequency  Domain  Framework  for  Analysis  and  Control  of 
Dynamical  Systems  in  Presence  of  Disturbances 

Andrzej  Banaszuk,  Pratt  &  Whitney  Program  Office,  UTRC 

9:00  am:  Newton  Iterations  for  Persistent  Periodic  Oscillations  in  Dynamical  Systems 
(KAM  Theory),  Umesh  Vaidya,  UCSB 

9:30  am:  Break 

Design  Of  Robust  Dynamical  Systems 

10:00  am:  Symmetry  vs  Robustness  of  Dynamical  Systems,  Prashant  Mehta,  UTRC 

10:30  am:  Aeroengine  Combustion  Stabilization,  Gregory  Hagen,  UTRC 

11:00  pm:  Coordinated  Control  of  UAVs,  Troy  Smith,  Shawn  Shadden,  California  Institute  of 
Technology 

11:30  pm:  DNA  example,  Igor  Mezic,  Professor,  UCSB,  Igor  Fedichenia,  UTRC 

12:00:  Working  Lunch 

1:00  pm:  Coarse  Methods  Analysis  of  Complex  Systems 

Yannis  Kevrekidis,  Professor,  Dept.  Chemical  Engineering,  Princeton  University 


